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{1+F^{-1}k^{2}}\]

with \(\alpha = -1, 0, 1\) and where F is the deformation radius. 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 & -1 & F^{-1}\partial_{x} \\ 1 & 0 & 0 \\ F^{-1}\partial_{x} & 0 & 0 \end{array}\right]\end{split}\]

which becomes, in Fourier space:

\[\begin{split}L = \left[\begin{array}{ccc} 0 & -1 & 2\pi i F^{-1} k_{1}/L \\ 1 & 0 & 0 \\ 2\pi i F_{-1} k_{1}/L & 0 & 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(control)[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.

Attributes

  • 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
  • _Froude : Froude number
  • _eps : Epsilon, used in QG formulation
  • _F : Rossby deformation radius

Methods

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

Parameters

  • control : Control object containing the relevant parameters/constants for initialisation
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,:,:]
create_eigenbasis(k)[source]

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

Parameters

-k1, k2 : The wavenumber pair for which the eigenvalue is desired

Returns

-A : The computed eigenvector array

create_eigenvalues(k)[source]

Creates the eigenvalues required for the creation of the eigenbasis (of the linear operator). Note that these have been found analytically and so an arbitrary L is not currently supported.

The eigenvalues returned are those for the QG formulation of the RSWE.

Parameters

-k1, k2 : The wavenumber pair for which the eigenvalue is desired

Returns

-eVals : The computed eigenvalues