User guide

This guide explains how the qgs framework can be used and configured to obtain the various model version.

1. Rationale behind qgs

The purpose of qgs is to provide a Python callable representing the model’s tendencies to the user, e.g. to provide a function \(\boldsymbol{f}\) (and also its Jacobian matrix) that can be integrated over time:

\[\dot{\boldsymbol{x}} = \boldsymbol{f}(t, \boldsymbol{x})\]

to obtain the model’s trajectory. This callable can of course also be used to perform any other kind of model analysis.

This guide provides:

  1. The different ways to configure qgs in order to obtain the function \(\boldsymbol{f}\) are explained in the section 2. Configuration of qgs.

  2. Some specific ways to configure qgs detailed in the section 3. Using user-defined symbolic basis.

  3. Examples of usages of the model’s tendencies function \(\boldsymbol{f}\) are given in the section 4. Using qgs (once configured).

2. Configuration of qgs

The computational flow to compute the function \(\boldsymbol{f}\) of the model’s tendencies \(\dot{\boldsymbol{x}} = \boldsymbol{f}(t, \boldsymbol{x})\) is explained in Computational flow and sketched below:

../_images/compuflow.png

Sketch of the computational flow.

This flow can be implemented step by step by the user, or can be automatically performed by the functions create_tendencies(). This function takes a QgParams parameters object and return the function \(\boldsymbol{f}\) and its Jacobian matrix. Optionally, it can also return the byproducts of the tendencies generation process, i.e. the objects containing the inner products between the model’s spatial modes and the tensors of the model’s tendencies terms.

This section of the user guide explains how to configure the QgParams object to obtain the desired model version and model’s parameters.

2.1 Two different computational methods and three different geophysical components

Two different methods are available to configure qgs, related to two different methods of computing the inner products between its spatial basis functions:

  1. Symbolic: It is the main way to specify the basis. The user can specify an arbitrary basis of functions on which the model’s partial differential equations will be projected. Depending on the dimension of the model version and on the computing power available, it might take a while to initialize.

  2. Analytic: The old method to specify the model version, in which the basis of functions are composed only of Fourier modes described in Coupled ocean-atmosphere model (MAOOAM) and is thus limited in term of modelling capacity. It is nevertheless the fastest way since it relies on analytic formulas derived in [USER-DCDV16] to compute the basis functions.

There are also three different geophysical components presently available in qgs:

  1. Atmosphere: This component is mandatory and provides a two-layers quasi-geostrophic atmosphere. See Atmospheric component for more details.

  2. Ocean: This component provides a shallow-water active oceanic layer. See Oceanic component for more details.

  3. Ground: This component provides a simple model for the ground (orography + heat exchange). It is present by default if only the atmospheric component is defined, but then with only the orography being activated (no heat exchange).

The components needed by the user and their parameters have to be defined by instantiating a QgParams object. How to create this object, initialize these components and set the parameters of the model is the subject of the next sections.

2.2 Initializing qgs

The model initialization first step requires the creation of a QgParams object:

from qgs.params.params import QgParams
model_parameters = QgParams()

This object contains basically all the information needed by qgs to construct the inner products and the tendencies tensor of the model, which are in turn needed to produces finally the model’s function \(\boldsymbol{f}\).

The different components required by the user need then to be specified, by providing information about the basis of functions used to project the partial differential equations of qgs. As said before, two methods are available:

2.2.1 The symbolic method

With this method, the user has to provide directly the basis of functions of each component. These functions have to be symbolic function expressions, and should be provided using Sympy. This has to be done using a SymbolicBasis object, which is basically a list of Sympy functions.

The user can construct his own basis (see below) or use the various built-in Fourier basis provided with qgs: ChannelFourierBasis or BasinFourierBasis. In the latter case, convenient constructor functions have been defined to help the user get the Fourier basis: contiguous_basin_basis() and contiguous_channel_basis(). These functions create contiguous Fourier basis for two different kind of boundary conditions (a zonal channel or a closed basin) shown on the first figure in Coupled ocean-atmosphere model (MAOOAM).

Note

A contiguous Fourier basis means here that the Fourier modes are all present in the model up to a given maximum wavenumber in each direction (zonal and meridional). Hence one has only to specify the maximum wavenumbers (and the model’s domain aspect ratio) to these constructor functions. One can also create non-contiguous Fourier basis by specifying wavenumbers explicitly at the ChannelFourierBasis or BasinFourierBasis instantiation (see the section 3.1 A simple example for an example).

Once constructed, the basis has to be provided to the QgParams object by using dedicated methods: set_atmospheric_modes(), set_oceanic_modes() and set_ground_modes(). With the constructor functions, one can activate the mandatory atmospheric layer by typing

from qgs.basis.fourier import contiguous_channel_basis
basis = contiguous_channel_basis(2, 2, 1.5)
model_parameters.set_atmospheric_modes(basis)

where we have defined a channel Fourier basis up to wavenumber 2 in both directions and an aspect ratio of \(1.5\).

Note

Please note that the aspect ratio of the basis object provided to qgs is not very important, because it is superseded by the aspect ratio sets in the QgParams object.

To activate the ocean or the ground components, the user has simply to use the method set_oceanic_modes() and set_ground_modes(). Note that providing a oceanic basis of functions automatically deactivate the ground component, and vice-versa.

Finally, since the MAOOAM Fourier basis are used frequently in qgs, convenient methods of the QgParams object allow one to create easily these basis inside this object (without the need to create them externally and then pass them to the qgs parameters object). These are the methods set_atmospheric_channel_fourier_modes(), set_oceanic_basin_fourier_modes() and set_ground_channel_fourier_modes(). For instance, the effect obtained with the 3 previous lines of code (activating the atmosphere) can also be obtained by typing:

model_parameters.set_atmospheric_channel_fourier_modes(2, 2, mode='symbolic')

These convenient methods can also initialize qgs with another method (called analytic) and which is described in the next section.

Warning

If you initialize one component with the symbolic method, all the other component that you define must be initialized with the same method.

2.2.2 The analytic method

Computing the inner products of the symbolic functions defined with Sympy can be very resources consuming, therefore if the basis of functions that you intend to use are the ones described in Coupled ocean-atmosphere model (MAOOAM), you might be interested to use the analytic method, which uses the analytic formula for the inner products given in [USER-DCDV16]. This initialization mode is put in action by using the convenient methods of the QgParams object: set_atmospheric_channel_fourier_modes(), set_oceanic_basin_fourier_modes() and set_ground_channel_fourier_modes().

For instance, to initialize a channel atmosphere with up to wavenumber 2 in both directions, one can simply write:

model_parameters.set_atmospheric_channel_fourier_modes(2, 2, mode='analytic')

Note that it is the default mode, so removing the mode argument will result in the same behavior.

Warning

If you initialize one component with the analytic method, all the other component that you define must be initialized with the same method.

2.3 Changing the default parameters of qgs

Now, how to change the parameters of qgs? As stated in the The model’s parameters module section of the References, there are seven types of parameters arranged in classes:

These parameters classes are regrouped into the global structure QgParams and are accessible through the attributes:

The parameters inside these structures can be changed by passing a dictionary of the new values to the set_params() method. For example, if one wants to change the Coriolis parameter \(f_0\) and the static stability of the atmosphere \(\sigma\), one has to write:

model_parameters.set_params({'f0': 1.195e-4, 'sigma':0.14916})

where model_parameters is an instance of the QgParams class. This method will find where the parameters are stored and will perform the substitution. However, some parameters may not have a unique name, for instance there is a parameter T0 for the stationary solution \(T_0\) of the 0-th order temperature for both the atmosphere and the ocean. In this case, one need to find out which part of the structure the parameter belongs to, and then call the set_params() of the corresponding object. For example, changing the parameter T0 in the atmosphere can be done with:

model_parameters.atemperature_params.set_params({'T0': 280.})

Finally, some specific methods allow to setup expansion 1 coefficients. Presently these are the AtmosphericTemperatureParams.set_thetas, AtmosphericTemperatureParams.set_insolation, OceanicTemperatureParams.set_insolation, GroundTemperatureParams.set_insolation and GroundParams.set_orography methods. For example, to activate the Newtonian cooling, one has to write:

model_parameters.atemperature_params.set_thetas(0.1, 0)

which indicates that the first component 2 of the radiative equilibrium mean temperature should be equal to \(0.1\).

Note

Using both the atmospheric Newtonian cooling coefficients with AtmosphericTemperatureParams.set_thetas and the heat exchange scheme AtmosphericTemperatureParams.set_insolation together doesn’t make so much sense. Using the Newtonian cooling scheme is useful when one wants to use the atmospheric model alone, while using the heat exchange scheme is useful when the atmosphere is connected to another component lying beneath it (ocean or ground).

Similarly, one activates the orography by typing:

model_parameters.ground_params.set_orography(0.2, 1)

We refer the reader to the description of these methods for more details (just click on the link above to get there).

Once your model is configured, you can review the list of parameters by calling the method QgParams.print_params():

model_parameters.print_params()

2.4 Creating the tendencies function

Once you have configured your QgParams instance, it is very simple to obtain the model’s tendencies \(\boldsymbol{f}\) and the its Jacobian matrix \(\boldsymbol{\mathrm{Df}}\). Just pass it to the function create_tendencies():

from qgs.functions.tendencies import create_tendencies

f, Df = create_tendencies(model_parameters)

The function f() hence produced can be used to generate the model’s trajectories. See the section 4. Using qgs (once configured) for the possible usages.

2.5 Saving your model

The simplest way to save your model is to pickle the functions generating the model’s tendencies and the Jacobian matrix. Hence, using the same name as in the previous section, one can type:

import pickle

# saving the model
model={'f': f, 'Df': Df, 'parameters': model_parameters}

with open('model.pickle', "wb") as file:
    pickle.dump(model, file, pickle.HIGHEST_PROTOCOL)

and it can be loaded again by typing

from qgs.params.params import QgParams

# loading the model
with open('model.pickle', "rb") as file:
    model = pickle.load(file)

f = model['f']
model_parameters = model['parameters']

Warning

Due to several different possible reasons, loading models saved previously on another machine may not work. The only thing to do is then to recompute the model tendencies with the loaded model parameters (using the function create_tendencies(). In this case, it is better to save only the model parameters:

import pickle

# saving the model

with open('model_parameters.pickle', "wb") as file:
    pickle.dump(model_parameters, file, pickle.HIGHEST_PROTOCOL)

It is also possible to save the inner products and/or the tensor storing the terms of the model’s tendencies. For instance, the function create_tendencies() allows to obtain these information:

f, Df, inner_products, tensor = create_tendencies(model_parameters, return_inner_products=True, return_qgtensor=True)

The objects QgParams, the inner products, and the object QgsTensor hence obtained can be saved using pickle or the built-in save_to_file() methods (respectively QgParams.save_to_file(), AtmosphericInnerProducts.save_to_file() and QgsTensor.save_to_file()).

Using these objects, it is possible to reconstruct by hand the model’s tendencies (see the section 3.2 A more involved example: Manually setting the basis and the inner products definition for an example).

3. Using user-defined symbolic basis

3.1 A simple example

In this simple example, we are going to create an atmosphere-ocean coupled model as in [USER-VDDCG15], but with some atmospheric modes missing.

First, we create the parameters object:

from qgs.params.params import QgParams
model_parameters = QgParams({'n': 1.5})

and we create a ChannelFourierBasis with all the modes up to wavenumber 2 in both directions, except the one with wavenumbers 1 and 2 in respectively the \(x\) and \(y\) direction:

from qgs.basis.fourier import ChannelFourierBasis
atm_basis = ChannelFourierBasis(np.array([[1,1],
                                          [2,1],
                                          [2,2]]),1.5)
model_parameters.set_atmospheric_modes(atm_basis)

Finally, we add the same ocean version as in [USER-VDDCG15]:

model_parameters.set_oceanic_basin_fourier_modes(2,4,mode='symbolic')

The last step is to set the parameters according to your needs (as seen in the section 2.3 Changing the default parameters of qgs).

The model hence configured can be passed to the function creating the model’s tendencies \(\boldsymbol{f}\), as detailed in the section 2.4 Creating the tendencies function.

3.2 A more involved example: Manually setting the basis and the inner products definition

Warning

This initialization method is not yet well-defined in qgs. It builds the model block by block to construct an ad-hoc model version.

In this section, we describe how to setup a user-defined basis for one of the model’s component. We will do it for the ocean, but the approach is similar for the other components. We will project the ocean equations on four modes proposed in [USER-Pie11]:

\[\begin{split}\tilde\phi_1(x,y) & = & 2\, e^{-\alpha x} \, \sin(\frac{n}{2} x)\, \sin(y), \nonumber \\ \tilde\phi_2(x,y) & = & 2\, e^{-\alpha x} \, \sin(n x)\, \sin(y), \nonumber \\ \tilde\phi_3(x,y) & = & 2\, e^{-\alpha x} \, \sin(\frac{n}{2} x)\, \sin(2 y), \nonumber \\ \tilde\phi_4(x,y) & = & 2\, e^{-\alpha x} \, \sin(n x)\, \sin(2 y), \nonumber \\\end{split}\]

and connect it to the channel atmosphere defined in the sections above, using a symbolic basis of functions. First we create the parameters object and the atmosphere:

from qgs.params.params import QgParams
model_parameters = QgParams({'n': 1.5})
model_parameters.set_atmospheric_channel_fourier_modes(2, 2, mode="symbolic")

Then we create a SymbolicBasis object:

from qgs.basis.base import SymbolicBasis
ocean_basis = SymbolicBasis()

We must then specify the function of the basis using Sympy:

from sympy import symbols, sin, exp
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))

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

ocean_basis.substitutions.append(('al', 1.))

The basis of function hence defined needs to be passed to the model’s parameter object:

model_parameters.set_oceanic_modes(ocean_basis)

and the user can set the parameters according to it needs (as seen in the previous section).

Additionally, for these particular basis function, a special inner product needs to be defined instead of the standard one proposed in Coupled ocean-atmosphere model (MAOOAM). We consider thus as in [USER-Pie11] and [USER-VDeCruz14] the following inner product:

\[(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\]

The special inner product can be specified in qgs by creating a user-defined subclass of the SymbolicInnerProductDefinition class defining the expression of the inner product:

from sympy import pi
from qgs.inner_products.definition import StandardSymbolicInnerProductDefinition
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)

and passing it to a OceanicSymbolicInnerProducts object:

ip = UserInnerProductDefinition()

from qgs.inner_products.symbolic import AtmosphericSymbolicInnerProducts, OceanicSymbolicInnerProducts

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

It will compute the inner products and may take a certain time (depending on your number of cpus available). Once computed, the corresponding tendencies must then be created manually, first by creating the QgsTensor object:

from qgs.tensors.qgtensor import QgsTensor

aotensor = QgsTensor(model_parameters, aip, oip)

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

from numba import njit
from qgs.functions.sparse_mul import sparse_mul3
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:]

This concludes the initialization of qgs, the function f() hence produced can be used to generate the model’s trajectories. An example of the construction exposed here, along with plots of the trajectories generated, can be found in the section Manual setting of the basis. See also the following section for the possible usages.

4. Using qgs (once configured)

Once the function \(\boldsymbol{f}\) giving the model’s tendencies has been obtained, it is possible to use it with the qgs built-in integrator to obtain the model’s trajectories:

from qgs.integrators.integrator import RungeKuttaIntegrator
import numpy as np

integrator = RungeKuttaIntegrator()
integrator.set_func(f)

ic = np.random.rand(model_parameters.ndim)*0.1  # create a vector of initial conditions with the same dimension as the model
integrator.integrate(0., 3000000., 0.1, ic=ic, write_steps=1)
time, trajectory = integrator.get_trajectories()

Note that it is also possible to use other ordinary differential equations integrators available on the market, see for instance the Example of DiffEqPy usage.

4.1 Analyzing the model’s output with diagnostics

In addition to the plotting of the model’s variables (the coefficients of the spectral expansion of the fields), it is also possible to plot several diagnostics associated to these variables, whether they are fields or scalars. These diagnostics are arranged in classes (Diagnostic being the parent abstract base class) and must be instantiated to obtain objects that convert the model output to the corresponding diagnostics.

For instance, we can define a MiddleAtmosphericStreamfunctionDiagnostic diagnostic that returns the 500hPa streamfunction field:

psi_a = MiddleAtmosphericStreamfunctionDiagnostic(model_parameters)

and feed it with a model outputs time and trajectory previously computed (like in the previous section):

psi_a.set_data(time, trajectory)

It is then possible to plot the corresponding fields with the plot() method:

psi_a.plot(50)
../_images/maooam_run-Copy1_63_1.png

or to easily create animation widget and movie of the diagnostic with respectively the animate() and movie() methods:

psi_a.movie(figsize=(13,8), anim_kwargs={'interval': 100, 'frames':100})

All the diagnostics classes should in principle have these 3 methods (plot, animate and movie). Note that the movies can also be saved in mpeg 4 format.

One can also define multiple diagnostics and couple them toghether with a MultiDiagnostic object. For instance, here we couple the previous psi_a diagnostic with a dynamical scatter plot:

variable_nondim = VariablesDiagnostic([2, 1, 0], model_parameters, False)
m = MultiDiagnostic(1,2)
m.add_diagnostic(variable_nondim, diagnostic_kwargs={'show_time':True, 'style':'2Dscatter'}, plot_kwargs={'ms': 0.2})
m.add_diagnostic(psi_a, diagnostic_kwargs={'show_time':False})
m.set_data(time, trajectory)

As with the other diagnostics, the MultiDiagnostic object can plot and create movies:

m.plot(50)
../_images/maooam_run-Copy1_70_0.png 
m.movie(figsize=(15,8), anim_kwargs={'interval': 100, 'frames':100})

Further examples using the diagnostics can be found in the sections Recovering the result of Reinhold and Pierrehumbert (1982), Recovering the result of De Cruz, Demaeyer and Vannitsem (2016) and Manual setting of the basis. Examples are also provided in the in the notebooks folder.

Moreover, the diagnostics classes methods are fully described in the section Diagnostics of the References.

Full list of currently available diagnostics

More diagnostics will be implemented soon.

4.2 Toolbox

The toolbox regroups submodules to make specific analysis with the model and are available in the qgs.toolbox module. For the references of these submodules, see Toolbox module.

Presently, the list of submodules available is the following:

  • lyapunov: A module to compute Lyapunov exponents and vectors.

More submodules will be implemented soon.

4.2.1 Lyapunov toolbox

This module allows you to integrate the model and simultaneously obtain the local Lyapunov vectors and exponents. The methods included are:

  1. The Benettin algorithm to compute the Forward and Backward Lyapunov vectors [USER-BGGS80]. This method is implemented in the LyapunovsEstimator class.

  2. The method of the intersection between the subspaces spanned by the Forward and Backward vectors to find the Covariant Lyapunov vectors [USER-ER85] (see also [USER-DPV21]). Suitable for long trajectories. This method is implemented in the CovariantLyapunovsEstimator class.

  3. The Ginelli method [USER-GPT+07] to compute also the Covariant Lyapunov vectors. Suitable for a trajectory not too long (depends on the memory available). This method is also implemented in the CovariantLyapunovsEstimator class.

See also [USER-KP12] for a description of these methods.

Some example notebooks on how to use this module are available in the notebooks/lyapunov folder.

5. Developers information

5.1 Running the test

The model core tensors can be tested by running pytest in the main folder:

pytest

This will run all the tests and return a report. The test cases are written using unittest. Additionally, test cases can be executed separately by running:

python -m unittest model_test/test_name.py

E.g., testing the MAOOAM inner products can be done by running:

python -m unittest model_test/test_inner_products.py

References

USER-BGGS80

Giancarlo Benettin, Luigi Galgani, Antonio Giorgilli, and Jean-Marie Strelcyn. Lyapunov characteristic exponents for smooth dynamical systems and for hamiltonian systems; a method for computing all of them. part 1: theory. Meccanica, 15(1):9–20, 1980. URL: https://link.springer.com/article/10.1007/BF02128236.

USER-DCDV16(1,2)

L. De Cruz, J. Demaeyer, and S. Vannitsem. The Modular Arbitrary-Order Ocean-Atmosphere Model: MAOOAM v1.0. Geoscientific Model Development, 9(8):2793–2808, 2016. URL: https://www.geosci-model-dev.net/9/2793/2016/.

USER-DPV21

Jonathan Demaeyer, Stephen G Penny, and Stéphane Vannitsem. Identifying efficient ensemble perturbations for initializing subseasonal-to-seasonal prediction. arXiv preprint arXiv:2109.07979, 2021. URL: https://arxiv.org/abs/2109.07979.

USER-ER85

J-P Eckmann and David Ruelle. Ergodic theory of chaos and strange attractors. The theory of chaotic attractors, pages 273–312, 1985. URL: https://link.springer.com/chapter/10.1007/978-0-387-21830-4_17.

USER-GPT+07

Francesco Ginelli, Pietro Poggi, Alessio Turchi, Hugues Chaté, Roberto Livi, and Antonio Politi. Characterizing dynamics with covariant lyapunov vectors. Physical review letters, 99(13):130601, 2007. URL: https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.99.130601.

USER-KP12

Pavel V Kuptsov and Ulrich Parlitz. Theory and computation of covariant lyapunov vectors. Journal of nonlinear science, 22(5):727–762, 2012. URL: https://link.springer.com/article/10.1007/s00332-012-9126-5.

USER-Pie11(1,2)

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. URL: https://journals.ametsoc.org/doi/full/10.1175/JPO-D-10-05018.1.

USER-VDeCruz14

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. URL: https://gmd.copernicus.org/articles/7/649/2014/.

USER-VDDCG15(1,2)

S. Vannitsem, J. Demaeyer, L. De Cruz, and M. Ghil. Low-frequency variability and heat transport in a low-order nonlinear coupled ocean–atmosphere model. Physica D: Nonlinear Phenomena, 309:71–85, 2015. URL: https://www.sciencedirect.com/science/article/abs/pii/S0167278915001335.

Footnotes

1

Generally in term of the specified basis functions.

2

The component corresponding to the first basis function of the atmopshere.