Exponential Integrator

This module implements the solution to the exponential integrator for the linear operator of the rotating shallow water equations as used by pypint.

Consider the RSWE, neglecting the nonlinear terms. We may write this as:

\[\frac{\partial u}{\partial t} = Lu\]

We may solve this via an integrating factor method, yielding:

\[u_{t} = e^{tL}u_{0}\]

to which we may freely apply various choices of timestepping. The exponential integrator then requires the efficient computation of the matrix exponential. We write the matrix exponential in the form:

\[e^{tL} = r_{k}^{\alpha} e^{t\omega_{k}^{\alpha}} (r_{k}^{\alpha})^{-1}\]

where \(r_{k}^{\alpha}\) is a matrix containing the eigenvectors of the linear operator (for a corresponding wavenumber, k, as we are working in Fourier space) and \(\omega_{k}^{\alpha}\) is a vector containing the associated eigenvalues of the linear operator, such that:

\[\omega_{k}^{\alpha} = \alpha\sqrt{f_{0}^{2} + 2\pi g H_{0} |k| / L}\]

with \(\alpha = -1, 0, 1\). This reduces the problem to an easier problem of finding the eigenvalues/eigenvectors and applying these. These are pre-computed for speed as they will not change over the course of the computation for constant gravity, rotation, and water depth.

Classes

-ExponentialIntegrator : Implements the exponential integrator for the RSWE

Notes

We write the linear operator, L, as:

\[\begin{split}L = \left[\begin{array}{ccc} 0 & -f_{0} & g\partial_{x} \\ f_{0} & 0 & g\partial_{y} \\ H_{0}\partial_{x} & H_{0}\partial_{y} & 0 \end{array}\right]\end{split}\]

which becomes, in Fourier space:

\[\begin{split}L = \left[\begin{array}{ccc} 0 & -f_{0} & 2\pi i g k_{1}/L \\ f_{0} & 0 & 2\pi i g k_{2}/L \\ 2\pi i H_{0} k_{1}/L & 2\pi i H_{0} k_{2}/L & 0 \end{array}\right]\end{split}\]

See Also

numpy

Authors: Terry Haut, Adam G. Peddle
Contact: ap553@exeter.ac.uk
Version: 2.0
class rswe_exponential_integrator.ExponentialIntegrator[source]

Implements the exponential integrator for the Rotating Shallow Water Equations.

This class implements the exponential integrator objects, which precompute and store the eigenvalues and eigenvectors of the linear operator in Fourier space and apply the matrix exponential for a given timestep when invoked. See module header for explanation of computation.

Methods

  • call : Invoke the matrix exponential
  • project_basis : Use the computed eigenbasis to project
call(U_hat, t)[source]

Call to invoke the exponential integrator on a given initial solution over some time.

Propagates the solution to the linear problem for some time, t, based on the initial condition, U. U is, in general, the value at the beginning of a timestep and t is the timestep, although this need not be the case.

Parameters

  • U_hat : the Fourier modes of the unknown solution, ordered as below
  • t : the time at which the solution is desired (or timestep to be used)

Returns

  • U_hat_sp : The computed solution, propagated by time t

Notes

  1. Works in -t direction as well.
  2. Input order: v1[:,:] = U_hat[0,:,:], v2[:,:] = U_hat[1,:,:], h[:,:] = U_hat[2,:,:]
project_basis(U_hat, t)[source]

Call to invoke the exponential integrator on a given initial solution over some time.

Propagates the solution to the linear problem for some time, t, based on the initial condition, U. U is, in general, the value at the beginning of a timestep and t is the timestep, although this need not be the case.

Parameters

  • U_hat : the Fourier modes of the unknown solution, ordered as below
  • t : the time at which the solution is desired (or timestep to be used)

Returns

  • U_hat_sp : The computed solution, propagated by time t

Notes

  1. Works in -t direction as well.
  2. Input order: v1[:,:] = U_hat[0,:,:], v2[:,:] = U_hat[1,:,:], h[:,:] = U_hat[2,:,:]
class rswe_exponential_integrator.ExponentialIntegrator_FullEqs(control)[source]

Implements the exponential integrator for the dimensional, perturbation-height Rotating Shallow Water Equations which have been symmetrised with respect to the geopotential height.

This class implements the exponential integrator objects, which precompute and store the eigenvalues and eigenvectors of the linear operator in Fourier space and apply the matrix exponential for a given timestep when invoked. See module header for explanation of computation.

Attributes

  • f0 : the Coriolis parameter
  • g : the gravitational acceleration
  • H0 : the mean water depth
  • deriv_factor : factor of 2*pi/L due to spectral differentiation
  • _Nx : Number of grid points, such that domain is Nx X Nx
  • eVals : Vector of eigenvalues, ordered alpha = 0, -1, 1
  • eBasis : Corresponding orthonormalised matrix of eigenvectors, arranged in columns

Methods

  • create_eigenvalues : Set up the analytically-computed eigenvalues for L
  • create_eigenbasis : Set up the analyticall-computed eigenbasis for L
  • call : Invoke the matrix exponential

Parameters

  • control : Control object containing the relevant parameters/constants for initialisation
create_eigenbasis(k1, k2)[source]

Create the eigenbasis of the linear operator. As with the eigenvalues, it is an implementation of an analytical solution.

This method computes the eigenvector matrix for a given wavenumber pair. It does not solve the eigenvalue problem, rather the analytical solution has been found for the representation of the linear operator which has been used.

Parameters

  • k1, k2 : the wavenumbers in the x- and y- directions

Returns

  • A : eigenvectors of L in columns, with slow mode 1st then fast waves in - and + directions

Notes

  1. Note that these have been found analytically and so an arbitrary L is not currently supported.
  2. The eigenbasis is orthonormalised
create_eigenvalues(k1, k2)[source]

Creates the eigenvalues required for the creation of the eigenbasis (of the linear operator).

This method returns the eigenvectors of the linear operator of the RSWE as formulated in spectral space. These have been found analytically to be:

\[evals = \alpha\sqrt{f_{0}^{2} + \beta^{2}gH_{0}|k|^{2}}\]

where \(\beta = 2\pi/L\) and \(\alpha = 0, -1, +1\), respectively.

Parameters

  • k1, k2 : the wavenumbers in the x- and y- directions

Returns

  • eVals : eigenvalues for slow mode and fast waves in - and + directions, in that order

Notes

Note that these have been found analytically and so an arbitrary L is not currently supported.