# Oceanic Basin Modes: Quasi-Geostrophic approach # =============================================== # # .. rst-class:: emphasis # # This tutorial was contributed by Christine Kaufhold and Francis # Poulin __. # # As a continuation of the Quasi-Geostrophic (QG) model described in the # other tutorial, we will now see how we can use Firedrake to compute # the spatial structure and frequencies of the freely evolving modes in this system, # what are referred to as basin modes. # Oceanic basin modes are low frequency structures that propagate # zonally in the oceans that alter the dynamics of Western Boundary Currents, # such as the Gulf Stream. In this particular tutorial we will show how to # solve the QG eigenvalue problem with no basic state and no dissipative # forces. # Unlike the other demo that integrated the equations forward in time, in this # problem it is necessary to compute the eigenvalues and eigenfunctions # for a particular differential operator. This requires using # PETSc __ matrices # and eigenvalue solvers in SLEPc __. # # This demo requires SLEPc and slepc4py to be installed. This is most easily # achieved by providing the optional --slepc flag to either firedrake-install # (for a new installation), or firedrake-update (to add SLEPc to an existing # installation). # # # Governing PDE # ~~~~~~~~~~~~~ # # We first briefly recap the nonlinear, one-layer QG equation that we # :doc:considered previously . # The interested reader can find the # derivations in :cite:QGeval-Pedlosky:1992 and :cite:QGeval-Vallis:2006. # This model consists of an evolution equation # for the Potential Vorticity, :math:q, and an elliptic problem through # which we can determine the streamfunction, # # .. math:: # # \partial_{t}q + \vec{\nabla}\cdot (\vec{u}q) + \beta v = 0 \\ # q = \nabla^{2} \psi - F\psi # # Where :math:\psi is the stream-function, :math:\vec{u}=(u, v) is the # velocity field, :math:q is the Potential Vorticity (PV), :math:\beta is the # Coriolis parameter and :math:F is the rotational Froude number. The velocity # field is easily obtained using # # .. math:: # # \vec{u} = \vec{\nabla}^{\bot}\psi, # \quad \mbox{ with } \quad # \vec{\nabla}^{\bot} = \hat{e_{z}} \times \vec{\nabla} # # We assume that the amplitude of the wave motion is very small, which # allows us to linearize the equations of motion and therefore neglect the # nonlinear advection, # # .. math:: \frac{\partial}{\partial t} (\nabla^{2} \psi - F\psi) = - \beta \frac{\partial \psi}{\partial x} # # We look for wave-like solutions that are periodic in time, with a # frequency of :math:\omega # # .. math:: \psi = \hat{\psi}(x, y)e^{-i\omega t} # # This has the advantage of removing the time derivative from the equation # and replacing it with an eigenvalue, :math:i \omega. By substituting # the above solution into the QG equation, we can find a complex # eigenvalue problem of the form # # .. math:: i\omega (\nabla^{2} \hat{\psi} - F\hat{\psi}) = \hat{\beta} \frac{\partial \hat{\psi}}{\partial x} # # Weak Formulation # ---------------- # # To use a finite element method it is necessary to formulate the weak # form and then we can use SLEPc in Firedrake to compute eigenvalue # problems easily. # To begin, we multiply this equation by a Test Function :math:\phi # and integrate over the domain :math:A. # # .. math:: # # i\omega \iint_{A} \Big(\phi\cdot\nabla^{2} \hat{\psi}\,dA - F\phi\hat{\psi}\,dA\Big) = \hat{\beta}\iint_{A} \phi \cdot \frac{\partial \hat{\psi}}{\partial x}\,dA # # To remove the Laplacian operator we use integration by parts and the Divergence theorem to obtain # # .. math:: # # \iint_{A} \phi \cdot \nabla^{2}\hat{\psi} \,dA = - \iint_{A} \nabla\phi \cdot \nabla\hat{\psi}\,dA + \oint_{\partial A} \phi \cdot \frac{\partial \hat{\psi}}{\partial n} \,dS # # No-normal flow boundary conditions are required and mathematically this # means that the streamfunction must be a constant on the boundary. Since # the test functions inherit these boundary conditions, # :math:\hat{\phi} = 0 on the boundary, the boundary integral # vanishes and the weak form becomes, # # .. math:: # # i\omega \iint_{A} \Big( \nabla\phi\cdot\nabla \hat{\psi}\,dA + F\phi\hat{\psi}\Big)\,dA = \hat{\beta}\iint_{A} \phi \cdot \frac{\partial \hat{\psi}}{\partial x}\,dA # # Firedrake code # -------------- # # Using this form, we can now implement this eigenvalue problem in # Firedrake. We import the Firedrake, PETSc, and SLEPc libraries. :: from firedrake import * from firedrake.petsc import PETSc try: from slepc4py import SLEPc except ImportError: import sys warning("Unable to import SLEPc, eigenvalue computation not possible (try firedrake-update --slepc)") sys.exit(0) # We specify the geometry to be a square geometry with :math:50 cells # with length :math:1. :: Lx = 1. Ly = 1. n0 = 50 mesh = RectangleMesh(n0, n0, Lx, Ly, reorder=None) # Next we define the function spaces within which our solution will # reside. :: Vcg = FunctionSpace(mesh,'CG',3) # We impose zero Dirichlet boundary conditions, in a strong sense, which # guarantee that we have no-normal flow at the boundary walls. :: bc = DirichletBC(Vcg, 0.0, "on_boundary") # The two non-dimensional parameters are the :math:\beta parameter, set # by the sphericity of the Earth, and the Froude number, the relative # importance of rotation to stratification. :: beta = Constant('1.0') F = Constant('1.0') # Additionally, we can create some Functions to store the eigenmodes. :: eigenmodes_real, eigenmodes_imag = Function(Vcg), Function(Vcg) # We define the Test Function :math:\phi and the Trial Function # :math:\psi in our function space. :: phi, psi = TestFunction(Vcg), TrialFunction(Vcg) # To build the weak formulation of our equation we need to build two PETSc # matrices in the form of a generalized eigenvalue problem, # :math:A\psi = \lambda M\psi. We impose the boundary conditions on the # mass matrix :math:M, since that is where we used integration by parts. :: a = beta*phi*psi.dx(0)*dx m = -inner(grad(psi), grad(phi))*dx - F*psi*phi*dx petsc_a = assemble(a).M.handle petsc_m = assemble(m, bcs=bc).M.handle # We can declare how many eigenpairs, eigenfunctions and eigenvalues, we # want to find :: num_eigenvalues = 1 # Next we will impose parameters onto our eigenvalue solver. The first is # specifying that we have an generalized eigenvalue problem that is # nonhermitian. The second specifies the spectral transform shift factor # to be non-zero. The third requires we use a Krylov-Schur method, # which is the default so this is not strictly necessary. Then, we ask for # the eigenvalues with the largest imaginary part. Finally, we specify the # tolerance. :: opts = PETSc.Options() opts.setValue("eps_gen_non_hermitian", None) opts.setValue("st_pc_factor_shift_type", "NONZERO") opts.setValue("eps_type", "krylovschur") opts.setValue("eps_largest_imaginary", None) opts.setValue("eps_tol", 1e-10) # Finally, we build our eigenvalue solver using SLEPc. We add our PETSc # matrices into the solver as operators and use setFromOptions() to call # the PETSc parameters we previously declared. :: es = SLEPc.EPS().create(comm=COMM_WORLD) es.setDimensions(num_eigenvalues) es.setOperators(petsc_a, petsc_m) es.setFromOptions() es.solve() # Additionally we can find the number of converged eigenvalues. :: nconv = es.getConverged() # We now get the real and imaginary parts of the eigenvalue and # eigenvector for the leading eigenpair (that with the largest in # magnitude imaginary part). First we check if we actually managed to # converge any eigenvalues at all. :: if nconv == 0: import sys warning("Did not converge any eigenvalues") sys.exit(0) # If we did, we go ahead and extract them from the SLEPc eigenvalue # solver:: vr, vi = petsc_a.getVecs() lam = es.getEigenpair(0, vr, vi) # and we gather the final eigenfunctions :: eigenmodes_real.vector()[:], eigenmodes_imag.vector()[:] = vr, vi # We can now list and show plots for the eigenvalues and eigenfunctions # that were found. :: print("Leading eigenvalue is:", lam) try: import matplotlib.pyplot as plt fig, axes = plt.subplots() colors = tripcolor(eigenmodes_real, axes=axes) fig.colorbar(colors) fig, axes = plt.subplots() colors = tripcolor(eigenmodes_imag, axes=axes) fig.colorbar(colors) except ImportError: warning("Matplotlib not available, not plotting eigemodes") # Below is a plot of the spatial structure of the real part of one of the eigenmodes computed above. # # .. figure:: eigenmode_real.png # :align: center # # Below is a plot of the spatial structure of the imaginary part of one of the eigenmodes computed above. # # .. figure:: eigenmode_imag.png # :align: center # # This demo can be found as a Python script in qgbasinmodes.py __. # # .. rubric:: References # # .. bibliography:: demo_references.bib # :filter: docname in docnames # :keyprefix: QGeval-