Introduction to FEM & Firedrake

Sia Ghelichkhan

Back to gadopt.org

- / -

Logo

Finite Element Method

A 15-Minute Introduction for G-ADOPT workshops

Sia Ghelichkhan and the G-ADOPT team

Objectives

By the end of this workshop, we want to have a (very) basic understanding of the finite element method:

  • Grasp the core concept: How continuous problems become discrete systems
  • Follow the mathematical journey: From PDEs to weak forms to matrix equations

Forward Modelling in Earth Sciences

Forward modelling is the process of using a known set of physical laws (mathematical equations), and input parameters to simulate or predict the behaviour of a physical system.

In essence, it involves solving a governing equation or system of equations to determine the state of the system under specified conditions.

Earth Systems
Any Earth system requiring forward modelling

Examples in Earth Sciences

Groundwater Groundwater Flow

Groundwater fills the spaces between soil particles and fractured rock beneath the earth’s surface.

Mantle Convection Mantle Convection

Mantle flow drives plate tectonics over geological timescales.

Glacial Isostasy Glacial Isostatic Adjustment

Loading and unloading of ice sheets causes crustal deformation and rebound.

Common challenge: Complex geometries, heterogeneous materials, and coupled physics require numerical methods like FEM.

Mathematical Formulation of Forward Problems

Forward modelling as a functional:

\mathbf{u} = \mathcal{F}(\mathbf{m})

Where:

  • \mathbf{u} are the state variables (temperature, velocity, pressure, displacement)
  • \mathbf{m} are the control parameters (material properties, boundary conditions, forcing terms)
  • \mathcal{F} is the forward operator (solves the governing PDE)

The functional relationship: \mathcal{F}: \mathbf{m} \mapsto \mathbf{u}

This maps control parameters to state variables through the solution of differential equations.

Examples:

Heat equation: \frac{\partial u}{\partial t} - \nabla \cdot (k \nabla u) = f

Stokes flow: -\nabla p + \nabla \cdot (\eta \nabla \mathbf{v}) = \mathbf{f}

Wave equation: \frac{\partial^2 u}{\partial t^2} - c^2 \nabla^2 u = f

Key insight: All involve derivatives → need numerical methods!

Core FEM Concepts

Mesh: Discretising Space

Mesh Progression
From continuous domain to discrete domain

Domain discretisation:

  1. Divide the domain into simple shapes (triangles, rectangles, cubes, …)
  2. Nodes are placed at vertices
  3. Elements are the individual triangles/tetrahedra

Key insight: The mesh quality affects solution accuracy. Finer mesh = better approximation.

The FEM Philosophy: “Just Good Enough!”

Strong form (exact solution): \nabla^2 u - f = 0 \quad \text{everywhere in } \Omega

A claim: u_{sol} is the solution \nabla^2 u_{sol} - f = 0

Problem: Too restrictive for complex domains!

Strong vs Weak

Test the claim with test functions v_i:

v_0(\nabla^2 u_{sol} - f) = 0v_1(\nabla^2 u_{sol} - f) = 0
v_2(\nabla^2 u_{sol} - f) = 0v_3(\nabla^2 u_{sol} - f) = 0

The complexity of our solution depends on the complexity of test functions. (which FunctionSpace we choose)

From Testing to Integration

The residual approach: \text{Res} = \int_\Omega v (\nabla^2 u_{sol} - f) \, dx = 0

Apply integration by parts: \int_\Omega \nabla v \cdot \nabla u_{sol} \, dx + \int_\Omega v f \, dx = 0

This becomes our discrete system: \mathbf{A} \mathbf{x} = \mathbf{b}

Why this works:

  • Relaxes differentiability requirements
  • Can handle discontinuous materials
  • Naturally leads to matrix systems
  • “Good enough” for engineering accuracy

Test and Trial Functions

Basis Functions
Hat functions: basis functions for linear elements

From infinite to finite dimensions:

Trial functions u_h: Our “claimed solution” u_h(x,y) = \sum_{i=1}^N U_i \phi_i(x,y)

Test functions v_h: The “testers” from our mesh

  • Choose v_h = \phi_j (one basis function at a time)
  • This gives us N equations for N unknowns U_i

Basis Functions

Basis Functions
Hat functions: basis functions for linear elements

Basis function properties:

  • Local support: \phi_i \neq 0 only near node i
  • Hat-shaped: Linear on each element
  • Interpolation: \phi_i(x_j) = \delta_{ij}

The connection:

  • Each mesh node → one basis function
  • Each basis function → one test function
  • Each test → one equation in \mathbf{A}\mathbf{x} = \mathbf{b}

Unified Form Language (UFL)

From Weak Form to Code

Remember our weak form: \int_\Omega \nabla v \cdot \nabla u_{sol} \, dx + \int_\Omega v f \, dx = 0

UFL makes this direct:

# Define function space on mesh
V = FunctionSpace(mesh, "Lagrange", 1)
u = TrialFunction(V)  # Our claimed solution
v = TestFunction(V)   # Our testers

# Define the weak form
a = dot(grad(v), grad(u))*dx  # Left side
L = v*f*dx                    # Right side

The magic of UFL:

  • Mesh → Function space
  • Basis functions → Trial/Test functions
  • Weak form → Bilinear/Linear forms
  • Assembly → Matrix system automatically!

2D Example: Heat Conduction in a Square

Problem Setup

2D Heat Problem
2D heat conduction in square domain

Physical problem: Heat conduction in a square

  • Domain: \Omega = [0,1] \times [0,1]
  • Governing equation: -\Delta u = f
  • Boundary conditions: u = 0 on boundary
  • Heat source: f = 10

UFL implementation:

# Mesh and function space
mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, "P", 1)

# Define problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(10.0)

# Variational forms
a = dot(grad(u), grad(v))*dx
L = f*v*dx