Idealised 2-D viscoelastic loading problem in an annulus¶

In this tutorial, we examine an idealised 2-D loading problem in an annulus domain.

This example focusses on differences between running simulations in a 2-D annulus and 2-D Cartesian domain. These can be summarised as follows:

1. The geometry of the problem - i.e. the computational mesh.
2. The radial direction of gravity (as opposed to the vertical direction in a Cartesian domain).
3. Solving a problem with laterally varying viscosity.
4. Accounting for a (rotational) velocity nullspace.

This example¶

Let's get started! The first step is to import the gadopt module, which provides access to Firedrake and associated functionality.

In [1]:

from gadopt import *
from gadopt.utility import step_func, vertical_component

from gadopt import * from gadopt.utility import step_func, vertical_component

We generate a circular manifold mesh (with 180 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 incremental displacement, we approximate the curved cylindrical shell domain quadratically, using the optional keyword argument degree$=2$.

As 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.

In [2]:

# Set up geometry:
rmin = 3480e3
rmax = 6371e3
D = rmax-rmin
nz = 32
ncells = 180
dz = D / nz

# Construct a circle mesh and then extrude into a cylinder:
surface_mesh = CircleManifoldMesh(ncells, radius=rmin, degree=2, name='surface_mesh')
mesh = ExtrudedMesh(surface_mesh, layers=nz, layer_height=dz, extrusion_type='radial')

bottom_id, top_id = "bottom", "top"
mesh.cartesian = False

# Set up geometry: rmin = 3480e3 rmax = 6371e3 D = rmax-rmin nz = 32 ncells = 180 dz = D / nz # Construct a circle mesh and then extrude into a cylinder: surface_mesh = CircleManifoldMesh(ncells, radius=rmin, degree=2, name='surface_mesh') mesh = ExtrudedMesh(surface_mesh, layers=nz, layer_height=dz, extrusion_type='radial') bottom_id, top_id = "bottom", "top" mesh.cartesian = False

We next set up the function spaces, and specify functions to hold our solutions.

In [3]:

# Set up function spaces - currently using the bilinear Q2Q1 element pair:
V = VectorFunctionSpace(mesh, "Q", 2)  # (Incremental) Displacement function space (vector)
W = FunctionSpace(mesh, "Q", 1)  # Pressure function space (scalar)
S = TensorFunctionSpace(mesh, "DQ", 2)  # (Discontinuous) Stress tensor function space (tensor)
R = FunctionSpace(mesh, "R", 0)  # Real function space (for constants)

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("Incremental Displacement")
z.subfunctions[1].rename("Pressure")

displacement = Function(V, name="displacement").assign(0)
stress_old = Function(S, name="stress_old").assign(0)

# Set up function spaces - currently using the bilinear Q2Q1 element pair: V = VectorFunctionSpace(mesh, "Q", 2) # (Incremental) Displacement function space (vector) W = FunctionSpace(mesh, "Q", 1) # Pressure function space (scalar) S = TensorFunctionSpace(mesh, "DQ", 2) # (Discontinuous) Stress tensor function space (tensor) R = FunctionSpace(mesh, "R", 0) # Real function space (for constants) 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("Incremental Displacement") z.subfunctions[1].rename("Pressure") displacement = Function(V, name="displacement").assign(0) stress_old = Function(S, name="stress_old").assign(0)

We can output function space information, for example the number of degrees of freedom (DOF).

In [4]:

# Output function space information:
log("Number of Incremental Displacement DOF:", V.dim())
log("Number of Pressure DOF:", W.dim())
log("Number of Velocity and Pressure DOF:", V.dim()+W.dim())

# Output function space information: log("Number of Incremental Displacement DOF:", V.dim()) log("Number of Pressure DOF:", W.dim()) log("Number of Velocity and Pressure DOF:", V.dim()+W.dim())

Number of Incremental Displacement DOF: 46800
Number of Pressure DOF: 5940
Number of Velocity and Pressure DOF: 52740

We can now visualise the resulting mesh.

In [5]:

active-ipynb

import pyvista as pv
import matplotlib.pyplot as plt

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)

import pyvista as pv import matplotlib.pyplot as plt 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)

Let's start initialising some parameters. First of all let's get the (symbolic) spatial coordinates of the mesh

In [6]:

X = SpatialCoordinate(mesh)

X = SpatialCoordinate(mesh)

Now we can set up the background profiles for the material properties. In this case the density and shear modulus vary in the vertical direction. We will approximate the series of layers using a smooth tanh function with a width of 40 km.

In [7]:

# layer properties from spada et al 2011
radius_values = [6371e3, 6301e3, 5951e3, 5701e3]
density_values = [3037, 3438, 3871, 4978]
shear_modulus_values = [0.50605e11, 0.70363e11, 1.05490e11, 2.28340e11]
viscosity_values = [2, -2, -2, -1.698970004]  # viscosity = 1e23 * 10**viscosity_values
# N.b. that we have modified the viscosity of the Lithosphere viscosity from
# Spada et al 2011 because we are using coarse grid resolution


def initialise_background_field(field, background_values, vertical_tanh_width=40e3):
    profile = background_values[0]
    sharpness = 1 / vertical_tanh_width
    depth = sqrt(X[0]**2 + X[1]**2)-radius_values[0]
    for i in range(1, len(background_values)):
        centre = radius_values[i] - radius_values[0]
        mag = background_values[i] - background_values[i-1]
        profile += step_func(depth, centre, mag, increasing=False, sharpness=sharpness)

    field.interpolate(profile)


density = Function(W, name="density")
initialise_background_field(density, density_values)

shear_modulus = Function(W, name="shear modulus")
initialise_background_field(shear_modulus, shear_modulus_values)

# layer properties from spada et al 2011 radius_values = [6371e3, 6301e3, 5951e3, 5701e3] density_values = [3037, 3438, 3871, 4978] shear_modulus_values = [0.50605e11, 0.70363e11, 1.05490e11, 2.28340e11] viscosity_values = [2, -2, -2, -1.698970004] # viscosity = 1e23 * 10**viscosity_values # N.b. that we have modified the viscosity of the Lithosphere viscosity from # Spada et al 2011 because we are using coarse grid resolution def initialise_background_field(field, background_values, vertical_tanh_width=40e3): profile = background_values[0] sharpness = 1 / vertical_tanh_width depth = sqrt(X[0]**2 + X[1]**2)-radius_values[0] for i in range(1, len(background_values)): centre = radius_values[i] - radius_values[0] mag = background_values[i] - background_values[i-1] profile += step_func(depth, centre, mag, increasing=False, sharpness=sharpness) field.interpolate(profile) density = Function(W, name="density") initialise_background_field(density, density_values) shear_modulus = Function(W, name="shear modulus") initialise_background_field(shear_modulus, shear_modulus_values)

Let's have a quick look at the density field using pyvista.

In [8]:

active-ipynb

# Read the PVD file
density_file = VTKFile('density.pvd').write(density)
reader = pv.get_reader("density.pvd")
data = reader.read()[0]  # MultiBlock mesh with only 1 block

# Create a plotter object
plotter = pv.Plotter(shape=(1, 1), border=False, notebook=True, off_screen=False)

# Make a colour map
boring_cmap = plt.get_cmap("viridis", 25)
# Add the warped displacement field to the frame
plotter.add_mesh(
    data,
    component=None,
    lighting=False,
    show_edges=False,
    cmap=boring_cmap,
    scalar_bar_args={
        "title": 'Density (kg / m^3)',
        "position_x": 0.8,
        "position_y": 0.2,
        "vertical": True,
        "title_font_size": 20,
        "label_font_size": 16,
        "fmt": "%.0f",
        "font_family": "arial",
    }
)
plotter.camera_position = 'xy'
plotter.show(jupyter_backend="static", interactive=False)
# Closes and finalizes movie
plotter.close()

# Read the PVD file density_file = VTKFile('density.pvd').write(density) reader = pv.get_reader("density.pvd") data = reader.read()[0] # MultiBlock mesh with only 1 block # Create a plotter object plotter = pv.Plotter(shape=(1, 1), border=False, notebook=True, off_screen=False) # Make a colour map boring_cmap = plt.get_cmap("viridis", 25) # Add the warped displacement field to the frame plotter.add_mesh( data, component=None, lighting=False, show_edges=False, cmap=boring_cmap, scalar_bar_args={ "title": 'Density (kg / m^3)', "position_x": 0.8, "position_y": 0.2, "vertical": True, "title_font_size": 20, "label_font_size": 16, "fmt": "%.0f", "font_family": "arial", } ) plotter.camera_position = 'xy' plotter.show(jupyter_backend="static", interactive=False) # Closes and finalizes movie plotter.close()