Ctraj has been applied to a new method of dynamical tracer reconstruction:
Download a pdf version of this page: description.pdf
DESCRIPTION OF CTRAJ
Contents
Basic algorithm
Tranformed coordinates
Contour advection
Semi-Lagrangian advection
Efficiency
Basic algorithm
The only thing that CTRAJ does is integrate the following pair of ordinary differential equations:
where x=(x, y) is position,
t is time
and v = (u, v) is a known velocity field.
All other codes are built up from this base.
The codes take a set of initial points (which we will call the initial "field")
and integrate them forward in time.
The velocity field is gridded in a Cartesian space and linearly interpoled.
The trajectories are integrated with a Runge Kutta fourth order
scheme with a constant time step.
Return to top of page
Transformed coordinates
The codes are general enough that any velocity field can
be integrated, so long as it is gridded in a Cartesian space.
Since the velocities are encapsulated in container classes,
adaptation to a non-gridded format or analytical computation
of velocities will be quite straightforward.
However the main purpose of these codes is to advect velocities
over the surface of the Earth (e.g. output from GCMs, assimilation data,
etc.) thus they must be transformed to a suitable projection.
We use an azimuthal equidistant projection limited to one hemisphere.
This has the dual advantage of eliminating the singularity at the
pole and having relatively little variation in the area of each grid cell.
The transformations from standard longitude-latitude
(constant diameter spherical polar) coodinates are as follows:
where r is defined as follows:
where R is the radius of the Earth
and h is the pole upon which the transformation is centred:
The resulting space has the following metric:
while the reverse transformations are given as:
The velocity field transforms as follows:
where u is meridional wind and v is zonal wind.
To run the simulation over the entire surface of the Earth
we use two fields, one for the Northern hemisphere and one for the Southern.
These fit together like two halves of an eggshell.
Since equatorial crossings are relatively rare,
we only check for them at intervals (typical when the field is output).
For this reason, the velocity fields overlap slightly, as seen in the figure.
The following transformations convert a coordinate to its equivalent
representation in the opposite hemisphere:
Where the subscript h labels the hemisphere.
Return to top of page
Contour advection
Consider a blob of dye stirred by a moving fluid.
To first order it can modelled by tracing its outlines
and allowing these to move with the fluid.
Contour advection is an adaptive method of modelling the
contours or isolines of a passive tracer.
It works by following the trajectories from points at intervals along
the contour.
To maintain the integrity of the curve,
new points are added or removed as needed.
We want the points to lie at roughly constant intervals
along the curve even as it is stretched in many directions
at once.
Rather than use intervals of constant distance, we use fraction of arc:
where 
is the difference in path
between two points and c is the radius of curvature and for a Cartesian
space is calculated as follows:
The path along the curve s is first approximated as a set of piece-wise
continous line segments.
The resulting parametric representation of the contour
is interpolated using a cubic spline.
A cubic spline returns a set of second-order
derivatives with which we can calculate the curvature:
New points are added between pairs of points that trace
out either too large an arc or too great a distance.
Likewise, points are removed if they are too close to their neighbour.
Return to top of page
Semi-Lagrangian advection
A passive tracer simulation aims to solve the following
differential equation, termed the advection equation:
where q is the tracer field in Lagrangian space,
x0 is the Lagrangian coordinate
and S is the source term.
While S can accomodate other sources,
such as chemical rates, here we concern ourselves with
the diffusion term:
where D is the diffusivity tensor.
In an Eulerian framework, this becomes:
Since, in the absence of diffusion, tracer values are constant
along the trajectories, we can build a tracer simulation
using only the trajectory integrator.
The idea is to define an Eulerian grid and
run back-trajectories from each grid point.
Values of the tracer for the current field are
interpolated from the previous field.
Because it doesn't matter in what order the trajectories
are run (they are independent of the tracer field)
and because the interpolation defines a series of coefficients
mapping a tracer field from one time step to the next,
we write these coefficients into a sparse matrix.
The sparse matrix will have sides the size of
the number of tracer grid points.
A mapping from a two-dimensional tracer field
to a one-dimensional vector must be defined.
Expressing this mathematically:
where q(ti) is the tracer field at the ith
timestep and Ai is a matrix.
A semi-Lagrangian tracer simulation is not subject to the CFL
criterion and is unconditionally stable.
The interpolation step provides a small amount of diffusion,
thus it can be almost eliminated by using an Eulerian time step
(that controls the interpolation between fields,
not the trajectory integration) that is as long
as the integration time.
Return to top of page
Efficiency
The advection codes spend the majority of their time
interpolating the velocity field.
We have mitigated this in the case of the tracer simulation
by writing coefficients to an array of sparse matrices
thus the tracer simulation needs only be run once over a given
period of time.
References
- D. G. Dritschel (1988). "Contour surgery: A topological reconnection scheme
for extended integrations using contour dynamics." Journal of
Computational Physics 77, 240-266.
- W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, B. P. (1992).
Numerical Recipes in C, 2nd Edition. Cambridge University Press,
Cambridge, UK.
Return to top of page
Next: Manual