Generalised stability analysis

Generalised stability analysis is a powerful tool for investigating the stability of physical systems. For an introduction to the mathematics, see the chapter in the mathematical background.

Computing the GST

Once the record of the simulation has been created, it is possible to perform generalised stability analysis with one line of code:

gst = compute_gst(initial, final, nsv)

The initial and final variables define the input and output of the propagator \(L\) respectively, while nsv is the number of singular vectors requested. By default, the mass matrices of the input and output spaces are used to define the norms on the input and output spaces. Optionally, the user can specify other matrices to change the norms used for either of these spaces. For example, to replicate the default behaviour, one could write

gst = compute_gst(initial, final, nsv, ic_norm=inner(u, v)*dx)

where u and v are instances of the TrialFunction and TestFunction classes on the appropriate input function space.

This one call will derive the tangent linear and adjoint systems, construct a matrix-free representation of the propagator, and use this inside a Krylov-Schur iteration to solve the GST singular value problem. This computation may take many iterations of the tangent linear and adjoint systems. The solution of the singular value problem is achieved with the SLEPc library.

Using the GST

Once the GST has been computed, it may be used as follows:

for i in range(gst.ncv):
  (sigma, u, v, residual) = gst.get_gst(i, return_vectors=True, return_residual=True)

The ncv member of the gst contains the number of converged singular vectors. This may be less than, equal to, or greater than the requested number of singular vectors.

By default, get_gst only returns the growth rate \(\sigma\) associated with the computed singular triplet. To fetch the singular vectors, pass return_vectors=True. To compute the residual of the eigenvalue computation, pass return_residual=True.

Example

A complete example of a generalised stability analysis of the tutorial example is presented below.

from __future__ import print_function
from dolfin import *
from dolfin_adjoint import *

n = 30
mesh = UnitSquareMesh(n, n)
V = VectorFunctionSpace(mesh, "CG", 2)

ic = project(Expression(("sin(2*pi*x[0])", "cos(2*pi*x[1])"), degree=2),  V)
nu = Constant(0.0001, name="nu")

def main(ic):
    u = ic.copy(deepcopy=True, name="Velocity")
    u_next = Function(V, name="VelocityNext")
    v = TestFunction(V)

    timestep = Constant(0.01)

    F = (inner((u_next - u)/timestep, v)
       + inner(grad(u_next)*u_next, v)
       + nu*inner(grad(u_next), grad(v)))*dx

    bc = DirichletBC(V, (0.0, 0.0), "on_boundary")

    t = 0.0
    end = 0.1
    while (t <= end):
        solve(F == 0, u_next, bc)
        u.assign(u_next)
        t += float(timestep)

    return u

if __name__ == "__main__":
    u = main(ic)

    gst = compute_gst("Velocity", "Velocity", nsv=3)
    for i in range(gst.ncv):
        sigma = gst.get_gst(i)
        print("Growth rate of vector %s: %s" % (i, sigma))

Download the gst code.

This prints the following output:

$ python tutorial11.py
...
Growth rate of vector 0: 4.09880352035
Growth rate of vector 1: 3.20037673764
Growth rate of vector 2: 3.07821571572
Growth rate of vector 3: 3.06242628866