2D tracer transport

This demo shows how the Firedrake DG advection equation demo can be implemented in Thetis.

The test case is the classic cosine-bell–cone–slotted-cylinder advection test case of [LeV96]. The domain is the unit square \(\Omega=[0,1]^2\) and the velocity corresponds to the solid body rotation \(\vec{u} = (0.5 - y, x - 0.5)\).

As usual, we start by importing Thetis.

from thetis import *

Define a 40-by-40 mesh of squares.

mesh2d = UnitSquareMesh(40, 40, quadrilateral=True)

We will solve a pure advection problem in non-conservative form, with no hydrodynamics. Therefore, bathymetry is not actually important. We set an arbitrary postive value, as this is required by Thetis to construct the solver object.

P1_2d = FunctionSpace(mesh2d, "CG", 1)
bathymetry2d = Function(P1_2d)
bathymetry2d.assign(1.0)
solver_obj = solver2d.FlowSolver2d(mesh2d, bathymetry2d)

To activate the tracer functionality of the 2D model, we need to add tracer fields. This is done using the ModelOptions2d.add_tracer_2d method.

The ‘label’ identifies the field inside Thetis. It should not contain spaces and typically ends with ‘_2d’ for 2D problems. The ‘name’ exists for users to identify the field and may contain spaces. It will appear in the colourbar of any vtk outputs. Finally, the ‘filename’ is used when storing outputs, so cannot contain spaces. The usual Thetis convention is to use CamelCase with a trailing ‘2d’.

Source terms and diffusivity coefficients should also be provided through this interface. We have a pure advection problem with no diffusivity or source terms. However, such terms can be specified by replacing the None values below.

options = solver_obj.options
labels = ['tracer_2d']
names = ['Depth averaged tracer']
filenames = ['Tracer2d']
options.fields_to_export = labels
for label, name, filename in zip(labels, names, filenames):
    options.add_tracer_2d(label, name, filename, source=None, diffusivity=None)

As mentioned above, we are only solving the tracer equation, which can be specified by setting tracer_only = True.

options.tracer_only = True

We will run for time \(2\pi\) – a full rotation – using a strong stability preserving third order Runge-Kutta method (SSPRK33). For consistency with the Firedrake demo, Thetis’ automatic timestep computation functionality is switched off and the simulation time is split into 600 steps, giving a timestep close to the CFL limit.

options.tracer_timestepper_type = 'SSPRK33'
options.timestep = pi/300.0
options.simulation_end_time = 2*pi
options.simulation_export_time = pi/15.0
options.tracer_timestepper_options.use_automatic_timestep = False

For consistency with the Firedrake demo, we do not use stabilization or slope limiters, both of which are used by default in Thetis. Slope limiters are used to obtain non-oscillatory solutions.

options.use_lax_friedrichs_tracer = False
options.use_limiter_for_tracers = False

The background tracer value is imposed as an upwind inflow condition. In general, this would be a Function, but here we just use a Constant value.

solver_obj.bnd_functions['tracer_2d'] = {'on_boundary': {'value': Constant(1.0)}}

The velocity field is set up using a simple analytic expression.

vP1_2d = VectorFunctionSpace(mesh2d, "CG", 1)
x, y = SpatialCoordinate(mesh2d)
uv_init = interpolate(as_vector([0.5 - y, x - 0.5]), vP1_2d)

Now, we set up the cosine-bell–cone–slotted-cylinder initial condition. The first four lines declare various parameters relating to the positions of these objects, while the analytic expressions appear in the last three lines. This code is simply copied from the Firedrake version of the demo.

bell_r0, bell_x0, bell_y0 = 0.15, 0.25, 0.5
cone_r0, cone_x0, cone_y0 = 0.15, 0.5, 0.25
cyl_r0, cyl_x0, cyl_y0 = 0.15, 0.5, 0.75
slot_left, slot_right, slot_top = 0.475, 0.525, 0.85

bell = 0.25*(1+cos(pi*min_value(sqrt(pow(x-bell_x0, 2) + pow(y-bell_y0, 2))/bell_r0, 1.0)))
cone = 1.0 - min_value(sqrt(pow(x-cone_x0, 2) + pow(y-cone_y0, 2))/cone_r0, 1.0)
slot_cyl = conditional(
    sqrt(pow(x-cyl_x0, 2) + pow(y-cyl_y0, 2)) < cyl_r0,
    conditional(And(And(x > slot_left, x < slot_right), y < slot_top), 0.0, 1.0),
    0.0
)

We then declare the inital condition, q_init, to be the sum of these fields. Furthermore, we add 1 to this, so that the initial field lies between 1 and 2, rather than between 0 and 1. This ensures that we can’t get away with neglecting the inflow boundary condition. We also save the initial state so that we can check the \(L^2\)-norm error at the end.

q_init = interpolate(1.0 + bell + cone + slot_cyl, P1_2d)
solver_obj.assign_initial_conditions(uv=uv_init, tracer_2d=q_init)

Now we are in a position to run the time loop.

solver_obj.iterate()

Finally, we display the normalised \(L^2\) error, by comparing to the initial condition.

q = solver_obj.fields.tracer_2d
L2_err = sqrt(assemble((q - q_init)*(q - q_init)*dx))
L2_init = sqrt(assemble(q_init*q_init*dx))
print_output(L2_err/L2_init)

This tutorial can be dowloaded as a Python script here.

References

[LeV96]

Randall J. LeVeque. High-Resolution Conservative Algorithms for Advection in Incompressible Flow. SIAM Journal on Numerical Analysis, 33(2):627–665, 1996. doi:10.1137/0733033.