# This magic makes plots appear
%matplotlib inline
import matplotlib.pyplot as plt
import os
os.environ["GADOPT_LOGLEVEL"] = "WARNING"
# Load Firedrake on Colab
try:
import firedrake
except ImportError:
!wget "https://github.com/g-adopt/tutorials/releases/latest/download/firedrake-install-real.sh" -O "/tmp/firedrake-install.sh" && bash "/tmp/firedrake-install.sh"
import firedrake
2-D Annulus Case¶
Weak formulation¶
As outlined in the previous tutorial, the weak formulation of the equations governing mantle dynamics are:
$$\int_\Omega (\nabla \vec{v})\colon \mu \left[ \nabla \vec{u} + \left( \nabla \vec{u} \right)^T\right] \ dx - \int_{\Omega} \left( \nabla \cdot \vec{v}\right)\ p \ dx - \int_{\Omega} Ra_0\ T\ \vec{v}\cdot\hat{k} \ dx = 0 \ \text{ for all } v\in V,$$
$$ \int_\Omega w \nabla \cdot \vec{u} \ dx\ \text{ for all } v\in V,$$
$$ \int_\Omega q\frac{\partial T}{\partial t} \ dx + \int_\Omega q \vec{u}\cdot \nabla T \ dx + \int_\Omega \left(\nabla q\right) \cdot \left(\kappa \nabla T\right) \ dx = 0 \text{ for all } q\in Q.$$
This example analyses mantle flow, subject to these same governing equations, 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 domain of different dimensions.
- With free-slip boundary conditions on both boundaries, this case incorporates a velocity nullspace, as well as a pressure nullspace.
from gadopt import *
import numpy as np
import pyvista as pv # Used for plotting vtk output
We first need a mesh. We generate a circular manifold mesh (with 32 elements) and extrude in the radial direction, using the optional keyword argument extrusion_type
, forming 8 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$.
# Set up geometry:
rmin, rmax, ncells, nlayers = 1.22, 2.22, 32, 8
mesh1d = CircleManifoldMesh(ncells, radius=rmin, degree=2)
mesh = ExtrudedMesh(mesh1d, layers=nlayers, extrusion_type='radial')
bottom_id, top_id = "bottom", "top"
n = FacetNormal(mesh) # used later for diagnostic purposes
domain_volume = assemble(1 * dx(domain=mesh)) # used later for diagnostic purposes
We can now visualise the resulting mesh:
V = FunctionSpace(mesh, "CG", 1)
File("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 now need the function spaces associated with this mesh, test functions and functions to hold our solutions. We follow the same definitions as our previous tutorial, using the bilinear Q2-Q1 element pair for velocity and pressure and a Q2 discretisation for temperature.
V = VectorFunctionSpace(mesh, "CG", degree=2) # Velocity function space (vector)
W = FunctionSpace(mesh, "CG", degree=1) # Pressure function space (scalar)
Q = FunctionSpace(mesh, "CG", degree=2) # Temperature function space (scalar)
Z = MixedFunctionSpace([V, W]) # Mixed function space
v, w = TestFunctions(Z)
q = TestFunction(Q)
z = Function(Z) # a field over the mixed function space Z.
u, p = split(z) # returns symbolic UFL expression for u and p
T = Function(Q, name="Temperature")
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.
# Set up temperature field and initialise:
X = SpatialCoordinate(mesh)
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 0x778549b96500>, 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), 7)), 25)
We can now visualise this initial temperature field:
File("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)
The Rayleigh number for this problem is defined in addition to the initial timestep $\Delta t$. The viscosity and thermal diffusivity are left at their default values (both = 1). These constants are used to create an Approximation representing the physical setup of the problem (options include Boussinesq, Extended Boussinesq, Truncated Anelastic Liquid and Anelastic Liquid), and a Timestep Adaptor, for controlling the time-step length (via a CFL criterion) as the simulation advances in time.
Ra, delta_t = Constant(1e5), Constant(1e-6)
approximation = BoussinesqApproximation(Ra)
t_adapt = TimestepAdaptor(delta_t, u, V, maximum_timestep=0.1, increase_tolerance=1.5)
Boundary conditions for temperature are set to $T = 0$ at the surface ($r_{\text{max}}$) and $T = 1$ at the base ($r_{\text{min}}$). For velocity, we specify free‐slip conditions on both boundaries. We incorporate these weakly through the Nitsche approximation. This illustrates a key advantage of the G-ADOPT framework: the user only specifies that the normal component of velocity is zero and all required changes are handled under the hood.
temp_bcs = {
bottom_id: {"T": 1.0},
top_id: {"T": 0.0},
}
stokes_bcs = {
bottom_id: {"un": 0},
top_id: {"un": 0},
}
As noted above, with a free-slip boundary condition on both boundaries, one can add an arbitrary rotation of the form $(-y, x)=r\hat{\mathbf{\theta}}$ to the velocity solution (i.e. this case incorporates a velocity nullspace, as well as a pressure nullspace). These lead to null-modes (eigenvectors) for the linear system, rendering the resulting matrix singular. In preconditioned Krylov methods these null-modes must be subtracted from the approximate solution at every iteration. We do that below, setting up a nullspace object as we did in the previous tutorial, albeit speciying the rotational
keyword argument to be True. Once again, this removes the requirement for a user to configure these options, further simplifying the task of setting up a (valid) geodynamical simulation.
Given the increased computational expense (typically requiring more degrees of freedom) in a 2-D annulus domain, the G-ADOPT library defaults to iterative solver parameters (in 2-D Cartesian domains, the framework defaults to direct solver parameters). Our iterative solver setup is configured to use the GAMG preconditioner for the velocity block of the Stokes system, to which we must provide near-nullspace information, which consists of three rotational (x_rotV
, y_rotV
, z_rotV
) and three translational (nns_x
, nns_y
, nns_z
) modes (see Davies et al. GMD, 2022 for details).
Z_nullspace = create_stokes_nullspace(Z, closed=True, rotational=True)
Z_near_nullspace = create_stokes_nullspace(Z, closed=False, rotational=True, translations=[0, 1])
We next come to solving the variational problem, with solver objects for the energy and Stokes systems created. For the energy system we pass in the solution field T, velocity u, the physical approximation, time step, temporal discretisation approach
(i.e. implicit middle point, being equivalent to a Crank Nicolson scheme) and boundary conditions. For the Stokes system, we pass in the solution fields z, Temperature, the physical approximation, boundary condition, the nullspace and near nullspace objects. We also set the
optional cartesian
keyword argument to False, which ensures that the unit vector, $\hat{k}$, points radially, in the direction opposite to gravity. Solution of the two variational problems is undertaken by PETSc.
energy_solver = EnergySolver(T, u, approximation, delta_t, ImplicitMidpoint, bcs=temp_bcs)
stokes_solver = StokesSolver(z, T, approximation, bcs=stokes_bcs,
cartesian=False,
nullspace=Z_nullspace, transpose_nullspace=Z_nullspace,
near_nullspace=Z_near_nullspace)
We can now initiate the time-loop, with the Stokes and energy systems solved seperately. These solve
calls once again convert symbolic mathematics into computation. In the time loop, set here to run for 2000 time-steps, we compute the RMS velocity and surface Nusselt number (using definitions from Jarvis, 1993) for diagnostic purposes, and print these results every 50 timesteps.
no_timesteps = 2000
time = 0
for timestep in range(0, no_timesteps):
dt = t_adapt.update_timestep()
time += dt
stokes_solver.solve()
energy_solver.solve()
# Compute diagnostics:
u_rms = sqrt(assemble(dot(u, u) * dx)) * sqrt(1./domain_volume)
f_ratio = rmin/rmax
top_scaling = -1.3290170684486309 # log(f_ratio) / (1.- f_ratio)
bot_scaling = -0.7303607313096079 # (f_ratio * log(f_ratio)) / (1.- f_ratio)
nu_top = (assemble(dot(grad(T), n) * ds_t) / assemble(1 * ds_t(domain=mesh))) * top_scaling
nu_base = (assemble(dot(grad(T), n) * ds_b) / assemble(1 * ds_b(domain=mesh))) * bot_scaling
if timestep % 50 == 0:
print(timestep, dt, u_rms, nu_top, nu_base)
0 1.5e-06 38.92507225950837 1.3290063767567633 -0.7304021210334157
50 0.0005528165283473966 74.45695663082371 6.750082575416409 -5.626283297339365
100 0.0003429408575886271 125.40476359475024 8.527681074428118 -8.351619069752307
150 0.0003186014175215388 184.62247182343648 9.330317487386871 -8.726040175836465
200 0.0002966119848032323 202.19611286677682 9.597635601360258 -9.157254089891918
250 0.0002957382421802955 197.54236486131134 9.688557353018389 -9.321638209070928
300 0.00030158078361295753 196.67630002157446 9.716770117208819 -9.226855392692926
350 0.00029963469426320244 198.22176808096148 9.71707884995764 -9.244952361960403
400 0.00029980403843654684 197.852802063836 9.728938298646082 -9.24767811001328
450 0.0003007318522135277 197.80217891473768 9.738615403997358 -9.235302379901192
500 0.00030065462610868917 197.96795583989984 9.742374109958252 -9.23678529598332
550 0.0003007791981715992 197.95199397786945 9.74642505382959 -9.235504246236431
600 0.00030093975144502465 197.96873995242322 9.749531270641079 -9.233375800567293
650 0.0003009909295883626 197.99194859097267 9.751514687153106 -9.232869499786506
700 0.0003010502419357017 197.99732673496322 9.753109309098086 -9.232177201326135
750 0.0003010955477040915 198.00543851674806 9.754288859072545 -9.231604013968951
800 0.00030112489328227644 198.011471832158 9.755141955968211 -9.231266500238533
850 0.0003011481137046225 198.0149888751054 9.75578668615229 -9.230977905803837
900 0.00030116555086359063 198.0180680948003 9.756261758819281 -9.23076421414417
950 0.00030117746893488 198.02025850765943 9.756612375858216 -9.230610742088562
1000 0.00030118657888780266 198.02179397679836 9.756871794481915 -9.230491633998879
1050 0.00030119281959006534 198.0229543710012 9.75706169542327 -9.230402526115554
1100 0.00030119706994914615 198.02375960298485 9.757199165893965 -9.230334701922938
1150 0.00030119984475017404 198.02430201766828 9.757296726692168 -9.230281493657023
1200 0.0003012012573404289 198.02463677445832 9.757362925627946 -9.230238741519186
1250 0.0003012017730216664 198.02478784066156 9.757403729649479 -9.230202619610326
1300 0.00030120129177756393 198.02477396915106 9.757423111332974 -9.230169915430096
1350 0.0003012001178156825 198.02459434921562 9.757422845111916 -9.230137928047867
1400 0.00030119811268410236 198.02424087899956 9.757403673625479 -9.230104006672002
1450 0.00030119540425989214 198.0236859658148 9.757364394711251 -9.230065425205858
1500 0.0003011918885278161 198.0228952062697 9.757302904192779 -9.230019294509137
1550 0.00030118756105761255 198.02181349228223 9.757215373899797 -9.229962239433716
1600 0.00030118238143600327 198.0203739179492 9.757096799201358 -9.229890411243176
1650 0.00030117632929115907 198.01848712236722 9.756940522710227 -9.229799138776203
1700 0.00030116951884277826 198.01604367284284 9.756738270347117 -9.229682948656468
1750 0.0003011621899477745 198.0129086508894 9.756480056408106 -9.229535249685517
1800 0.00030115495941157 198.00891464799724 9.756153662410892 -9.229348253327682
1850 0.0003011489386272989 198.0038567611258 9.755744514896689 -9.229112667512565
1900 0.0003011461852572527 197.9974736979201 9.755234367637922 -9.228817350001819
1950 0.00030114156466427394 197.98942802632368 9.754599922276388 -9.228448633179426
We can now visualise the final temperature field:
File("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)