Brian E. J. Rose, University at Albany
This document uses the interactive Jupyter notebook format. The notes can be accessed in several different ways:
- The interactive notebooks are hosted on
githubat https://github.com/brian-rose/ClimateModeling_courseware - The latest versions can be viewed as static web pages rendered on nbviewer
- A complete snapshot of the notes as of May 2017 (end of spring semester) are available on Brian's website.
Also here is a legacy version from 2015.
Many of these notes make use of the climlab package, available at https://github.com/brian-rose/climlab
- The one-dimensional diffusion equation
- Discretizing the diffusion operator in space
- Coding the discretized diffusion operator in
numpy - Discretizing the time derivative
- Stability analysis of the FTCS scheme
- Numerical tests with a shorter timestep
- The need for a more efficient method
- Implicit time method
- Your homework assignment
Suppose that a quantity
The diffusive flux is
There will be local changes in
Putting this together gives the classical diffusion equation in one dimension
For simplicity, we are going to limit ourselves to Cartesian geometry rather than meridional diffusion on a sphere.
We will also assume here that
This equation represents a time-dependent diffusion process. It is an initial-boundary value problem. We want to integrate the model forward in time to model the changes in the field
Solving a differential equation on a computer always requires some approximation to represent the continuous function
We have already dealt with simple discretization of the time derivative back in [Lecture 2](Lecture02 -- Solving the zero-dimensional EBM.ipynb). We used the forward Euler method to step all our of radiation models forward in time so far.
We will discretize time and space on grids
so that
The governing equation can be written in terms of the convergence of the diffusive flux:
It is sensible to use a centered difference to approximate this derivative:
$$ \frac{\partial F}{\partial x} \bigg|j \approx \frac{F{j+\frac{1}{2}} - F_{j-\frac{1}{2}}}{x_{j+\frac{1}{2}} - x_{j-\frac{1}{2}}} $$
The time tendency at point
$$ \frac{\partial u}{\partial t} \bigg|j \approx - \frac{F{j+\frac{1}{2}} - F_{j-\frac{1}{2}}}{x_{j+\frac{1}{2}} - x_{j-\frac{1}{2}}} $$
The flux itself depends on a spatial derivative of
But we actually want to approximate
and
Putting this all together, we can write the time tendency at
$$ \frac{\partial u}{\partial t} \bigg|j \approx K \frac{ \frac{u{j+1} - u_{j}}{x_{j+1} - x_{j}} - \frac{u_{j} - u_{j-1}}{x_{j} - x_{j-1}}}{x_{j+\frac{1}{2}} - x_{j-\frac{1}{2}}} $$
We'll make things easy on ourselves by using uniform grid spacing in
So our final formula for the diffusive flux convergence is
$$ \frac{\partial u}{\partial t} \bigg|j \approx K \frac{ u{j+1} - 2 u_{j} + u_{j-1}}{\Delta x^2} $$
Suppose the domain is
The physical boundary condition at the walls is that there can be no flux in or out of the walls:
So the boundary conditions on
Suppose we have a grid of
$x^*_0 = 0 $ $x^*_1 = \Delta x$ $x^*_2 = 2~\Delta x$ - ...
$x^*_j = j~\Delta x$ - ...
$x^*_{J-1} = (J-1)~\Delta x = 1 - \Delta x $ $x^*_J = J ~ \Delta x = 1 $
Clearly then the grid spacing must be
We'll define the fluxes on this grid. The boundary conditions can thus be written
Since our centered difference discretization defines
The first grid point for
$x_0 = \Delta x / 2$ $x_1 = \Delta x / 2 + \Delta x$ $x_2 = \Delta x / 2 + 2~\Delta x$ - ...
$x_j = \Delta x / 2 + j~\Delta x$ - ...
$x_{J-1} = \Delta x / 2 + (J-1)~\Delta x = 1 - \Delta x / 2 $
At
Subbing in
The same procedure at the other wall yields
$$ \frac{\partial u}{\partial t} \bigg|{J-1} \approx - K \frac{ u{J-1} - u_{J-2} }{\Delta x^2} $$
Pulling this all together we have a complete discretization of the diffusion operator including its boundary conditions:
$$ \frac{\partial u}{\partial t} \bigg|j \approx K \frac{ u{j+1} - 2 u_{j} + u_{j-1}}{\Delta x^2}, ~~~~~~ j=1,...,J-2 $$
$$ \frac{\partial u}{\partial t} \bigg|{J-1} \approx - K \frac{ u{J-1} - u_{J-2} }{\Delta x^2} $$
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import display, Math, LatexHere we will divide our domain up into 20 grid points.
J1 = 20
J = J1
deltax = 1./J
display(Math(r'J = %i' %J))
display(Math(r'\Delta x = %0.3f' %deltax))The fluxes will be solved on the staggered grid with 21 points.
xstag = np.linspace(0., 1., J+1)
x = xstag[:-1] + deltax/2
print x[ 0.025 0.075 0.125 0.175 0.225 0.275 0.325 0.375 0.425 0.475
0.525 0.575 0.625 0.675 0.725 0.775 0.825 0.875 0.925 0.975]
u = np.zeros_like(x)Here's one way to implement the finite difference, using array indexing.
dudx = (u[1:] - u[:-1]) / (x[1:] - x[:-1])dudx.shape(19,)
We can also use the function numpy.diff() to accomplish the same thing:
help(np.diff)Help on function diff in module numpy.lib.function_base:
diff(a, n=1, axis=-1)
Calculate the n-th discrete difference along given axis.
The first difference is given by ``out[n] = a[n+1] - a[n]`` along
the given axis, higher differences are calculated by using `diff`
recursively.
Parameters
----------
a : array_like
Input array
n : int, optional
The number of times values are differenced.
axis : int, optional
The axis along which the difference is taken, default is the last axis.
Returns
-------
diff : ndarray
The n-th differences. The shape of the output is the same as `a`
except along `axis` where the dimension is smaller by `n`.
.
See Also
--------
gradient, ediff1d, cumsum
Examples
--------
>>> x = np.array([1, 2, 4, 7, 0])
>>> np.diff(x)
array([ 1, 2, 3, -7])
>>> np.diff(x, n=2)
array([ 1, 1, -10])
>>> x = np.array([[1, 3, 6, 10], [0, 5, 6, 8]])
>>> np.diff(x)
array([[2, 3, 4],
[5, 1, 2]])
>>> np.diff(x, axis=0)
array([[-1, 2, 0, -2]])
np.diff(u).shape(19,)
Here is a function that computes the diffusive flux
def diffusive_flux(u, deltax, K=1):
# Take the finite difference
F = np.diff(u)/deltax
# add a zero as the first element (no flux on boundary)
F = np.insert(F, 0, 0.)
# add another zero as the last element (no flux on boundary)
F = np.append(F, 0.)
# flux is DOWN gradient, proportional to D
return -K*Fdiffusive_flux(u,deltax).shape(21,)
The time tendency of
def diffusion(u, deltax, K=1):
# compute flux
F = diffusive_flux(u, deltax, K)
# take convergence of flux
return -np.diff(F) / deltaxSuppose we have an initial
The gaussian (bell curve) function is a convenient way to create such a field.
def gaussian(x, mean, std):
return np.exp(-(x-mean)**2/(2*std**2))/np.sqrt(2*np.pi*std**2)K = 0.01
u = gaussian(x, 0.5, 0.08)
dudt = diffusion(u, deltax, K=K)
fig, ax = plt.subplots(1)
ax.plot(x, u, label='$u(x)$')
ax.plot(x, dudt, label='$du/dt$')
ax.legend()<matplotlib.legend.Legend at 0x1067b4f50>
Hopefully this makes sense. The diffusion is acting to smooth out
Use a random number generator to create some noisy initial conditions.
fig = plt.figure(figsize=(10,8))
for n in range(4):
u = np.random.random(J)
dudt = diffusion(u, deltax, K)
ax = fig.add_subplot(2,2,n+1)
ax.plot(x, u)
ax.plot(x, dudt)
The simplest way to discretize the time derivative is the forward Euler method:
We have already used this method to step our prognostic variables forward in time.
Solving the above for the future value of
We apply our discretization of the diffusion operator to the current value of the field
(except at the boundaries, where the diffusion operator is slightly different -- see above).
Together, this scheme is known as Forward Time, Centered Space or FTCS.
It is very simple to implement in numpy code.
def step_forward(u, deltax, deltat, K=1):
dudt = diffusion(u, deltax, K)
return u + deltat * dudtK = 0.01
deltat = 0.125
deltat1 = deltat
u0 = gaussian(x, 0.5, 0.08)
u1 = step_forward(u0, deltax, deltat1, K)
fig, ax = plt.subplots(1)
ax.plot(x, u0, label='initial')
ax.plot(x, u1, label='next')
ax.legend()<matplotlib.legend.Legend at 0x11015afd0>
Let's loop through a number of timesteps.
# regular resolution
J = 20
deltax = 1./J
xstag = np.linspace(0., 1., J+1)
x = xstag[:-1] + deltax/2u = gaussian(x, 0.5, 0.08)
niter = 11
for n in range(niter):
u = step_forward(u, deltax, deltat1, K)
plt.plot(x, u, label=n)
plt.legend()<matplotlib.legend.Legend at 0x1101eda90>
The numerics were easy to implement, and the scheme seems to work very well! The results are physically sensible.
Try setting
What happens?
# double the resolution
scaling_factor = 2
J = J1 * scaling_factor
deltax = 1./J
xstag = np.linspace(0., 1., J+1)
x = xstag[:-1] + deltax/2u = gaussian(x, 0.5, 0.08)
for n in range(niter):
u = step_forward(u, deltax, deltat1, K)
plt.plot(x, u, label=n)
plt.legend()<matplotlib.legend.Legend at 0x10ff92fd0>
Suddenly our scheme is producing numerical noise that grows in time and overwhelms to smooth physical solution we are trying to model.
This is bad!
What went wrong, and what can we do about it?
Following Press et al. (1988), "Numerical Recipes in C: The Art of Scientific Computing", Cambridge University Press.
This is an example of the so-called von Neumann Stability Analysis. It is a form of normal mode analysis for a discrete system.
We look for normal mode solutions (i.e. wavy sines and cosines) of the finite difference equations of the form
where
The number
The question is, under what conditions do wavy solutions grow with time? (This is bad, as it means small numerical noise will become large numerical noise and make our differencing scheme unusable)
Let's substitute the normal mode solution into our finite difference equation
Divide through by $\xi^n \exp(ikj~\Delta x)$:
The exponentials simplify
Or using a double angle identity,
We need to prevent growing normal modes. So successive amplitudes should be
The stability condition is thus
and this condition must be met for EVERY possible wavenumber
Because
We conclude the the FTCS scheme is stable so long as this stability condition is met:
The maximum timestep we can use with the FTCS scheme for the diffusion equation is proportional to
** A doubling of the spatial resolution would require a 4x shorter timestep to preserve numerical stability. **
Physically, the restriction is that the maximum allowable timestep is approximately the diffusion time across a grid cell of width
Going back to our Gaussian example, let's double the resolution but shorten the timestep by a factor of 4.
# double the resolution
J = J1 * scaling_factor
deltax = 1./J
xstag = np.linspace(0., 1., J+1)
x = xstag[:-1] + deltax/2K = 0.01
# The maximum stable timestep
deltat_max = deltax**2 / 2 / K
print 'The maximum allowable timestep is %f' %deltat_max
deltat = deltat1 / scaling_factor**2
print '4x the previous timestep is %f' %deltatThe maximum allowable timestep is 0.031250
4x the previous timestep is 0.031250
u = gaussian(x, 0.5, 0.08)
for n in range(niter):
for t in range(scaling_factor**2):
u = step_forward(u, deltax, deltat, K)
plt.plot(x, u, label=n)
plt.legend()<matplotlib.legend.Legend at 0x110538110>
Success! The graph now looks like a smoother (higher resolution) version of our first integration with the coarser grid.
But at a big cost: our calculation required 4 times more timesteps to do the same integration.
The total increase in computational cost was actally a factor of 8 to get a factor of 2 increase in spatial resolution.
In practice the condition
is often too restrictive to be practical!
Consider our diffusive EBM. Suppose we want a spatial resolution of 1º latitude. Then we have 180 grid points from pole to pole, and our physical length scale is
We were using a diffusivity of $ D = 0.6 ~ \text{W m}^{-2}~\text{K}^{-1} $ and a heat capacity of
Accounting for the spherical geometry in our EBM, this translates to
Recall that this is the diffusivity associated with the large-scale motion of the atmosphere (mostly). If we take a typical velocity scale for a mid-latitude eddy,
Using these numbers the stability condition is roughly
which is less than one hour!
And if we wanted to double the resolution to 0.5º, we would need a timestep of just a few minutes.
This can be a very onerous requirement for a model that would like to integrate out for many years. We can do better, but we need a different time discretization!
With numerical methods for partial differential equations, it often turns out that a small change in the discretization can make an enormous difference in the results.
The implicit time scheme applies exactly the same centered difference scheme to the spatial derivatives in the diffusion operator.
But instead of applying the operator to the field
The scheme looks like
$$ \frac{u_j^{n+1} - u_j^n}{\Delta t} = \frac{K}{\Delta x^2} \left( u^{n+1}{j+1} - 2 u^{n+1}{j} + u^{n+1}_{j-1} \right) $$
in the interior, and at the boundaries:
and
This might seem like a strange way to write the system, since we don't know the future state of the system at
Let's move all terms evaluated at
$$ u_j^{n+1} - \frac{K \Delta t}{\Delta x^2} \left( u^{n+1}{j+1} - 2 u^{n+1}{j} + u^{n+1}_{j-1} \right) = u_j^n $$
or
$$ -K^* u^{n+1}{j+1} + \left(1+2K^* \right) u_j^{n+1} - K^* u{j-1}^{n+1} = u_j^n $$
(in the interior)
where we have introduced a non-dimensional diffusivity
We can write this as a matrix equation
where
and
Solving for the future state of the system
Solving a tridiagonal matrix problem like this is a very common operation in computer science, and efficient numerical routines are available in many languages (including Python / numpy!)
We'll skip the details, but the amplification factor for this scheme is (see Numerical Recipes book or other text on numerical methods):
so the stability criterion of $$ \bigg| \frac{\xi^{n+1}}{\xi^n} \bigg| \le 1 $$
is met for any value of
The implicit method (also called backward time) is unconditionally stable for any choice of timestep.
Write Python code to solve the diffusion equation using this implicit time method. Demonstrate that it is numerically stable for much larger timesteps than we were able to use with the forward-time method. One way to do this is to use a much higher spatial resolution.
We have just scratched the surface of the wonders and sorrows of numerical methods here. The implicit method is very stable but is not the most accurate method for a diffusion problem, particularly when you are interested in some of the faster dynamics of the system (as opposed to just getting the system quickly to its equilibrium state).
There are always trade-offs in the choice of a numerical method.
The equations for most climate models are sufficiently complex that more than one numerical method is necessary. Even in the simple diffusive EBM, the radiation terms are handled by a forward-time method while the diffusion term is solved implicitly.
Once you have worked through the above problem (diffusion only), you might want to look in the climlab code to see how the diffusion solver is implemented there, and how it is used when you integrate the EBM.
%load_ext version_information
%version_information numpy, matplotlib| Software | Version |
|---|---|
| Python | 2.7.12 64bit [GCC 4.2.1 (Based on Apple Inc. build 5658) (LLVM build 2336.11.00)] |
| IPython | 5.3.0 |
| OS | Darwin 16.5.0 x86_64 i386 64bit |
| numpy | 1.11.1 |
| matplotlib | 2.0.0 |
| Thu May 25 11:31:18 2017 EDT | |
The author of this notebook is Brian E. J. Rose, University at Albany.
It was developed in support of ATM 623: Climate Modeling, a graduate-level course in the Department of Atmospheric and Envionmental Sciences
Development of these notes and the climlab software is partially supported by the National Science Foundation under award AGS-1455071 to Brian Rose. Any opinions, findings, conclusions or recommendations expressed here are mine and do not necessarily reflect the views of the National Science Foundation.