Idealised 2-D mantle convection problem inside an annulus¶
In this tutorial, we analyse mantle flow in a 2-D annulus domain. We define our domain by the radii of the inner ($r_{\text{min}}$) and outer ($r_{\text{max}}$) boundaries. These are chosen such that the non-dimensional depth of the mantle, $z = r_{\text{max}} - r_{\text{min}} = 1$, and the ratio of the inner and outer radii, $f=r_{\text{min}} / r_{\text{max}} = 0.55$, thus approximating the ratio between the radii of Earth's surface and core-mantle-boundary (CMB). Specifically, we set $r_{\text{min}} = 1.22$ and $r_{\text{max}} = 2.22$.
This example focusses on differences between running simulations in a 2-D annulus and 2-D Cartesian domain. These can be summarised as follows:
- The geometry of the problem - i.e. the computational mesh.
- The radial direction of gravity (as opposed to the vertical direction in a Cartesian domain).
- Initialisation of the temperature field in a different domain.
- With free-slip boundary conditions on both boundaries, this case incorporates a (rotational) velocity nullspace, as well as a pressure nullspace.
The example is configured at $Ra = 1e5$. Boundary conditions are free-slip at the surface and base of the domain.
The first step is to import the gadopt module, which provides access to Firedrake and associated functionality. We also import pyvista, which is used for plotting vtk output.
from gadopt import *
import pyvista as pv
We next set up the mesh, function spaces, and specify functions to hold our solutions, as with our previous tutorials.
We generate a circular manifold mesh (with 128 elements) and extrude in the radial direction,
using the optional keyword argument extrusion_type, forming 32 layers. To better represent the
curvature of the domain and ensure accuracy of our quadratic representation of velocity, we
approximate the curved cylindrical shell domain quadratically, using the optional keyword argument degree$=2$.
Because this problem is not formulated in a Cartesian geometry, we set the mesh.cartesian
attribute to False. This ensures the correct configuration of a radially inward vertical direction.
rmin, rmax, ncells, nlayers = 1.22, 2.22, 128, 32
mesh1d = CircleManifoldMesh(ncells, radius=rmin, degree=2) # construct a circle mesh
mesh = ExtrudedMesh(mesh1d, layers=nlayers, extrusion_type='radial') # extrude into a cylinder
mesh.cartesian = False
bottom_id, top_id = "bottom", "top"
V = VectorFunctionSpace(mesh, "CG", 2) # Velocity function space (vector)
W = FunctionSpace(mesh, "CG", 1) # Pressure function space (scalar)
Q = FunctionSpace(mesh, "CG", 2) # Temperature function space (scalar)
Z = MixedFunctionSpace([V, W]) # Mixed function space.
z = Function(Z) # A field over the mixed function space Z.
u, p = split(z) # Returns symbolic UFL expression for u and p
z.subfunctions[0].rename("Velocity")
z.subfunctions[1].rename("Pressure")
We can now visualise the resulting mesh.
VTKFile("mesh.pvd").write(Function(V))
mesh_data = pv.read("mesh/mesh_0.vtu")
edges = mesh_data.extract_all_edges()
plotter = pv.Plotter(notebook=True)
plotter.add_mesh(edges, color="black")
plotter.camera_position = "xy"
plotter.show(jupyter_backend="static", interactive=False)
We next specify the important constants for this problem, and set up the approximation.
Ra = Constant(1e5) # Rayleigh number
approximation = BoussinesqApproximation(Ra)
As with the previous examples, we set up a Timestep Adaptor, for controlling the time-step length (via a CFL criterion) as the simulation advances in time. For the latter, we specify the initial time, initial timestep $\Delta t$, and number of timesteps. Given the low Rayleigh number, a steady-state tolerance is also specified, allowing the simulation to exit when a steady-state has been achieved.
time = 0.0 # Initial time
delta_t = Constant(1e-7) # Initial time-step
timesteps = 20000 # Maximum number of timesteps
t_adapt = TimestepAdaptor(delta_t, u, V, maximum_timestep=0.1, increase_tolerance=1.5)
steady_state_tolerance = 1e-7 # Used to determine if solution has reached a steady state.
We next set up and initialise our Temperature field. We choose the initial temperature distribution to trigger upwelling of 4 equidistant plumes. This initial temperature field is prescribed as:
$$T(x,y) = (r_{\text{max}} - r) + A\cos(4 \; atan2\ (y,x)) \sin(r-r_{\text{min}}) \pi)$$
where $A=0.02$ is the amplitude of the initial perturbation.
X = SpatialCoordinate(mesh)
T = Function(Q, name="Temperature")
r = sqrt(X[0]**2 + X[1]**2)
T.interpolate(rmax - r + 0.02*cos(4*atan2(X[1], X[0])) * sin((r - rmin) * pi))
Coefficient(WithGeometry(FunctionSpace(<firedrake.mesh.ExtrudedMeshTopology object at 0x73984eb0aba0>, TensorProductElement(FiniteElement('Lagrange', interval, 2), FiniteElement('Lagrange', interval, 2), cell=TensorProductCell(interval, interval)), name=None), Mesh(VectorElement(TensorProductElement(FiniteElement('Lagrange', interval, 2), FiniteElement('Lagrange', interval, 1), cell=TensorProductCell(interval, interval)), dim=2), 6)), 29)
We can plot this initial temperature field:
VTKFile("temp.pvd").write(T)
temp_data = pv.read("temp/temp_0.vtu")
plotter = pv.Plotter(notebook=True)
plotter.add_mesh(temp_data)
plotter.camera_position = "xy"
plotter.show(jupyter_backend="static", interactive=False)