Lecture 6 (b)

Unstructured grids

Most problems have complex domains.
Boundary conditions crucial.
Little “structure” to problem.
What use is “structure” in the grid?
Use unstructured grids.
But how, whilst keeping accuracy?

An unstructued grid representing parts of the UK

Function bases

A Fourier Series

\[ \phi = \sum_n c_n(t) \chi_n(x) = \sum_{n=-N}^N c_n(t) \exp(\mathrm{i} n x) \]

allow us to compute eg $\partial_x \phi$ everywhere.

Extend to more complex bases $\chi_n(x)$ to

match geometry better;
minimize Gibb’s oscillations;
speed up computations.

In the discussion of finite differencing methods, Fourier Series and Fourier analysis was mostly introduced as a way of studying stability and accuracy. Its use in spectral methods was mentioned, as was the idea of spectral accuracy, where the error decreases exponentially with the number of terms in the sum.

However, there are a range of issues with just using Fourier Series. First, they only have nice properties (like spectral accuracy) on periodic domains. This can make them useful on angles in spherical geometries, but not when boundaries are involved. Second, the well-known Gibb’s oscillations cause huge problems when the solution is discontinuous. Finally, the cost of performing the Fourier transforms (or similar) to compute the coefficients can be expensive.

Generally, most numerical schemes can be thought of as acting on a function basis of some sort, so it is worth spending some time looking at useful approaches.

Coefficients

How does finite $\{\chi_n \colon n=0,\dots,N\}$ best match solution to

\[ \partial_t \phi + F(\phi) = 0? \]

Case 1: match at points $x_j$ (colocation).

\[ c_n \chi_n(x_j)=\phi(x_j), \quad \forall x_j, \quad j=0,\dots,N. \]

Gives linear system for $c_n$.

Case 2: minimise residual $R(\phi)=\partial_t \phi+F(\phi)$ wrt basis (Galerkin).

\[ \int R(\phi) \chi_m(x) = 0, \quad \forall m = 0, \dots, N. \]

With a function basis the (approximate) solution is, in principle, known everywhere in space and time. That means we can compute the solution and its derivatives everywhere, and therefore check how well it matches the PDE of interest everywhere. It’s impossible to match it match perfectly, as the function basis is truncated to make it possible to calculate on a computer.

We have two choices. We either make it solve the PDE exactly at a finite number of points, or we make it as good as possible “on average”. The first approach is called colocation and has some significant advantages in understanding what is going on an computing terms when the PDE is highly nonlinear. However, we have no control on how the solution behaves between the colocation points, and the best choice of colocation points then becomes extremely important.

The Galerkin approach is equivalent to minimizing the error in the $L_2$ norm, and projects the residual with respect to the basis itself.

Galerkin methods

Galerkin methods best “on average”. PDE

\[ \partial_t \phi+F(\phi)=0 \]

becomes, for $\phi = \sum c_n(t) \chi_n(x)$,

\[ \begin{aligned} \sum_n I_{nk} \frac{\mathrm{d} c_n}{\mathrm{d} t} &= - \int F \left( \sum_n c_n \chi_n \right) \chi_k, \quad \forall k, \\ I_{nk} &= \int \chi_n \chi_k \, . \end{aligned} \]

Coupled nonlinear ODEs.

Spectral Convergence

Spectral convergence with $N$.
Faster than all $\Delta x^k \sim N^{-k}$ from eg finite differencing.
Timestep limited by “effective resolution”, so $\Delta t \sim N^{-2}$.
Problems with
- non-smooth solutions;
- non-smooth domains/boundaries.

A spectral method advects the initial data $\sin^4(2 \pi x)$ around the domain $x \in [0, 1]$ up to $t = 0.5$. Fourier modes $c_n = -N, \dots, N$ are used. Spectral (exponential) convergence with the number of modes is seen. It is also clear that this is much faster than first or second order convergence as indicated on the bottom panel.

Here we see the huge advantages of spectral methods. Even using only N=8 the solution is visually indistinguishable from the high resolution solution that we are treating as “exact”. That is not quite as small a number as you might think: the basis function expansion runs from -N to N so there are 2 N + 1 modes, and each mode is complex, rather than real. So N=8 corresponds to 34 floating point numbers.

The key advantage is the rate of convergence, which is exponential with N. This is faster than any polynomial convergence rate, and astonishingly quick compared to the first and second order convergence of standard finite difference methods, for example. This is sketched on the plot, where it is assumed that the equivalent first or second order method matches the accuracy of the spectral method when there are 4 degrees of freedom. We see that such methods are nowhere near as efficient.

There are a number of problems of varying importance. First, the CFL constraint is much stronger for spectral methods, going as the inverse square of the number of degrees of freedom N, rather than just the inverse. This is a serious problem if a large number of modes are needed. However, mostly they are not.

Much more significant are the problems that arise when either the solution or the domain is not a simple as this smooth, periodic case. The Gibbs’ effect was mentioned previously, and can really kick in here. Essentially, trying to represent a discontinuous function using a standard function basis as in a Fourier series does not give pointwise convergence near the discontinuity - there are oscillations there. This leads to a positive feedback loop when used in a PDE - these oscillations grow in an unstable fashion. A number of methods have been proposed to fix this - in my opinion none are satisfactory. Less severe problems occur when the solution is continuous but not smooth. This typically restricts the convergence rate to be polynomial, with the order linked to the degree of differentiability of the solution. Thus, for solutions that are rougher than $C^4$ it’s often the case that finite difference methods are more efficient, as the spectral methods will struggle to converge faster and are more computationally costly.

A wider problem is that spectral methods are extremely sensitive to the boundary conditions imposed, and the domain on which they are imposed. When the boundaries are periodic everything is fine; when the boundaries have small, rough features (such as an ocean or atmosphere model in the vertical direction) the spectral methods can have problems converging until the number of modes used is impractically high.

Summary

Semi-discretization simplifies constructing complex schemes.
Function basis methods give solution everywhere using truncated series approximation.
Can suffer Gibb’s effects at eg discontinuities.
Spectral methods minimize coupling of basis functions.
Spectral convergence exponential.
Spectral methods often too sensitive for practical use.