Quasi-Geostrophic Model

This tutorial was contributed by Francis Poulin, based on code from Colin Cotter.

The Quasi-Geostrophic (QG) model is very important in geophysical fluid dynamics as it describes some aspects of large-scale flows in the oceans and atmosphere very well. The interested reader can find derivations in [Ped92] and [Val06].

In these notes we present the nonlinear equations for the one-layer QG model with a free-surface. Then, the weak form will be derived as is needed for Firedrake.

Governing Equations

The Quasi-Geostrophic (QG) model is very similar to the 2D vorticity equation. Since the leading order geostrophic velocity is incompressible in the horizontal, the governing equations can be written as

\[\begin{split}\begin{aligned} \partial_t q + \vec \nabla \cdot \left( \vec u q \right) + \beta v &= 0, \\ \vec u & = \vec\nabla^\perp \psi, \\ \nabla^2 \psi - \frac{1}{L_d^2} \psi &= q. \end{aligned}\end{split}\]

where the \(\psi\) and \(q\) are the streamfunction and Potential Vorticity (PV). The Laplacian is 2D since we are only in the horizontal plane and we defined

\[\vec\nabla^\perp = \hat e_z \times \vec\nabla.\]

The first equation above states that the PV is conserved following the flow. The second equation forces the leading order velocity to be geostrophic and the third equation is the definition for the QG PV for this barotropic model. To solve this using Finite Elements it is necessary to establish the weak form of the model, which is done in the next subsection.

Weak Form

Evolving the nonlinear equations consists of two steps. First, the elliptic problem must be solved to compute the streamfunction given the PV. Second, the PV equation must be integrated forward in time. This is done using a strong stability preserving Runge Kutta 3 (SSPRK3) method.

Elliptic Equation

First, we focus on the elliptic inversion in the case of a flat bottom. If we compute the inner product of the equation with the test function \(\phi\) we obtain,

\[\begin{split}\begin{aligned} \langle \nabla^2 \psi, \phi \rangle - \frac{1}{L_d^2} \langle \psi, \phi \rangle &= \langle q, \phi \rangle, \\ \langle \nabla \psi, \nabla \phi \rangle + \frac{1}{L_d^2} \langle \psi, \phi \rangle &= -\langle q, \phi \rangle,\end{aligned}\end{split}\]

where in the second equation we used the divergence theorem and the homogeneous Dirichlet boundary conditions on the test function.

Evolution Equation

The SSPRK3 method used as explained in [Got05] can be written as

\[\begin{split}\begin{aligned} q^{(1)} &= q^n - \Delta t \left[ \vec \nabla \cdot \left( \vec u^n q^n \right) + \beta v^n \right] , \\ q^{(2)} &= \frac34 q^n + \frac14 \left[ q^{(1)} - \Delta t \vec \nabla \cdot \left( \vec u^{(1)} q^{(1)} \right) - \Delta t \beta v^{(1)}\right], \\ q^{n+1} &= \frac13 q^n + \frac23 \left[ q^{(2)} - \Delta t \vec \nabla \cdot \left( \vec u^{(2)} q^{(2)} \right) - \Delta t \beta v^{(1)} \right].\end{aligned}\end{split}\]

To get the weak form we need to introduce a test function, \(p\), and take the inner product of the first equation with \(p\).

\[\begin{split}\begin{aligned} \langle q^{(1)}, p \rangle &= \langle q^n, p \rangle - \Delta t \langle \vec \nabla \cdot \left( \vec u^n q^n \right), p \rangle - \Delta t \langle \beta v, q \rangle, \\ \langle q^{(1)}, p \rangle - \Delta t \langle \vec u^n q^n, \vec\nabla p \rangle + \Delta t \langle \beta v, q \rangle &= \langle q^n, p \rangle - \Delta t \langle \vec u^n q^n, p \rangle_{bdry}\end{aligned}\end{split}\]

The first and second terms on the left hand side are referred to as \(a_{mass}\) and \(a_{int}\) in the code. The first term on the right-hand side is referred to as \(a_{mass}\) in the code. The second term on the right-hand side is the extra term due to the DG framework, which does not exist in the CG version of the problem and it is referred to as \(a_{flux}\). This above problem must be solved for \(q^{(1)}\) and then \(q^{(2)}\) and then these are used to compute the numerical approximation to the PV at the new time \(q^{n+1}\).

We now move on to the implementation of the QG model for the case of a freely propagating Rossby wave. As ever, we begin by importing the Firedrake library.

from firedrake import *

Next we define the domain we will solve the equations on, square domain with 50 cells in each direction that is periodic along the x-axis.

Lx   = 2.*pi                                     # Zonal length
Ly   = 2.*pi                                     # Meridonal length
n0   = 50                                        # Spatial resolution
mesh = PeriodicRectangleMesh(n0, n0, Lx, Ly,  direction="x", quadrilateral=True)

We define function spaces:

Vdg = FunctionSpace(mesh,"DG",1)               # DG elements for Potential Vorticity (PV)
Vcg = FunctionSpace(mesh,"CG",1)               # CG elements for Streamfunction
Vu  = VectorFunctionSpace(mesh,"DG",1)          # DG elements for velocity

and initial conditions for the potential vorticity, here we use Firedrake’s ability to interpolate UFL expressions.

x = SpatialCoordinate(mesh)
q0 = Function(Vdg).interpolate(0.1*sin(x[0])*sin(x[1]))

We define some Functions to store the fields:

dq1 = Function(Vdg)       # PV fields for different time steps
qh  = Function(Vdg)
q1  = Function(Vdg)

psi0 = Function(Vcg)      # Streamfunctions for different time steps
psi1 = Function(Vcg)

along with the physical parameters of the model.

F    = Constant(1.0)         # Rotational Froude number
beta = Constant(0.1)         # beta plane coefficient
Dt   = 0.1                   # Time step
dt   = Constant(Dt)

Next, we define the variational problems. First the elliptic problem for the stream function.

psi = TrialFunction(Vcg)
phi = TestFunction(Vcg)

# Build the weak form for the inversion
Apsi = (inner(grad(psi),grad(phi)) + F*psi*phi)*dx
Lpsi = -q1*phi*dx

We impose homogeneous dirichlet boundary conditions on the stream function at the top and bottom of the domain.

bc1 = DirichletBC(Vcg, 0., (1, 2))

psi_problem = LinearVariationalProblem(Apsi,Lpsi,psi0,bcs=bc1)
psi_solver = LinearVariationalSolver(psi_problem,
                                     solver_parameters={
        'ksp_type':'cg',
        'pc_type':'sor'
        })

Next we’ll set up the advection equation, for which we need an operator \(\vec\nabla^\perp\), defined as a python anonymouus function:

gradperp = lambda u: as_vector((-u.dx(1), u.dx(0)))

For upwinding, we’ll need a representation of the normal to a facet, and a way of selecting the upwind side:

n = FacetNormal(mesh)
un = 0.5*(dot(gradperp(psi0), n) + abs(dot(gradperp(psi0), n)))

Now the variational problem for the advection equation itself.

q = TrialFunction(Vdg)
p = TestFunction(Vdg)
a_mass = p*q*dx
a_int  = (dot(grad(p), -gradperp(psi0)*q) + beta*p*psi0.dx(0))*dx
a_flux = (dot(jump(p), un('+')*q('+') - un('-')*q('-')) )*dS
arhs   = a_mass - dt*(a_int + a_flux)

q_problem = LinearVariationalProblem(a_mass, action(arhs,q1), dq1)

Since the operator is a mass matrix in a discontinuous space, it can be inverted exactly using an incomplete LU factorisation with zero fill.

q_solver  = LinearVariationalSolver(q_problem,
                                    solver_parameters={
        'ksp_type':'preonly',
        'pc_type':'bjacobi',
        'sub_pc_type': 'ilu'
        })

To visualise the output of the simulation, we create a File object. To which we can store multiple Functions. So that we can distinguish between them we will give them descriptive names.

q0.rename("Potential vorticity")
psi0.rename("Stream function")
v = Function(Vu, name="gradperp(stream function)")
v.project(gradperp(psi0))

output = File("output.pvd")

output.write(q0, psi0, v)

Now all that is left is to define the timestepping parameters and execute the time loop.

t = 0.
T = 10.
dumpfreq = 5
tdump = 0

v0 = Function(Vu)

while(t < (T-Dt/2)):

  # Compute the streamfunction for the known value of q0
  q1.assign(q0)
  psi_solver.solve()
  q_solver.solve()

  # Find intermediate solution q^(1)
  q1.assign(dq1)
  psi_solver.solve()
  q_solver.solve()

  # Find intermediate solution q^(2)
  q1.assign(0.75*q0 + 0.25*dq1)
  psi_solver.solve()
  q_solver.solve()

  # Find new solution q^(n+1)
  q0.assign(q0/3 + 2*dq1/3)

  # Store solutions to xml and pvd
  t += Dt
  print(t)

  tdump += 1
  if tdump == dumpfreq:
      tdump -= dumpfreq
      v.project(gradperp(psi0))
      output.write(q0, psi0, v, time=t)

A python script version of this demo can be found here.

References

[Got05]Sigal Gottlieb. On high order strong stability preserving Runge–Kutta and multi step time discretizations. Journal of Scientific Computing, 25(1):105–128, 2005. doi:10.1007/s10915-004-4635-5.
[Ped92]Joseph Pedlosky. Geophysical Fluid Dynamics. Springer study edition. Springer New York, 1992. ISBN 9780387963877.
[Val06]Geoffrey K. Vallis. Atmospheric and Oceanic Fluid Dynamics. Cambridge University Press, Cambridge, U.K., 2006.