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