sl.tracer.calculate.trajectories: Calculate kinematic trajectories

In FESOM/spheRlab: Spherical Geometry, Analysis, and Plotting of Geoscientific Data on Arbitrary Grids

Calculate kinematic trajectories

Description

A suite of three different integration methods for kinematic trajectory calculations on unstructured triangular meshes.

Usage

sl.tracer.calculate.trajectories(method = "Petterssen", U, x.ini, t.ini = 0, grid, i.neighs, tri.cont.ini = NULL, abort.if = NULL, tgr = NULL, dt.def = 1, T.end = 10, maxiter = 4, redfac = 2,  cart_geo= "geo", patch.level = 10, Rsphere = 6371)

Arguments

method	a string defining the integration scheme, either 'Petterssen', 'Euler' or 'Pansy'. The latter is currently only available for Cartesian coordinates.
U	a numerical array (dim: [2, N.nodes, N.time]), where N.time is the number of available time steps and N.nodes the number of nodes defining the grid. U[1,,] contains the velocity component in zonal (or x-) direction, and U[2,,] contains the velocity component in meridional (or y-) direction.
x.ini	a numerical array (dim: [2, N.ini]) containing the initial positions. Longitudes (or x-values) are therefore given by x.ini[1,] and latitudes (or y-values) by x.ini[2,]. N.ini is the number of initial positions.
t.ini	a numerical scalar, the time for which the initial positions x.ini are valid. Default is t.ini=0.
grid	a list containing at least the entries lon, lat and elem. The former two are numerical vectors containing the coordinates of the nodes (with length N.nodes), the latter is a numerical array (dim: [N.elem,3]) defining the FEM grid by storing the indices of the node coordinates in the rows.
i.neighs	a list of length N.elem which contains all elements (indices) that share at least one node with the given element, i.e. i.neighs[[17]] contains the indices of all neighboring elements for the element with index 17. If not provided, this list will be created during runtime.
tri.cont.ini	a numerical vector of length N.ini, containing the triangles (indices) that contain the initial positions. If not provided, it will be determined during runtime, which can be quite exhaustive for large numbers of elements.
abort.if	<documentation pending>
tgr	a numerical vector, the time grid of the input data. If not speficied, it is set to 0:(N.time-1). The other time-related variables must match the units of this vector.
dt.def	a numeric scalar defining the default time step size (in units of tgr).
T.end	a numeric scalar defining the maximal length of the time integration (in units of tgr).
maxiter	an integer scalar defining the maximum number of time step adjustments per time step
qfrac	a scalar number (greater than 1) by which the time step is divided during the adjustment
cart_geo	a string, either 'cart' or 'geo', defining if the grid shall be treated as geographic coordinates ('geo') or Cartesian coordinates ('cart').
geo.gc_rad	a string, either 'gc' or 'rad', defining how the displacement is calculated for geographical coordinates. Either as great-circle distance ('gc') or approximating with small angles in northward and eastward direction ('rad').
patch.level	an integer number defining a search radius for triangle tracking, I recommend keeping the default.
Rsphere	a numerical scalar, defining the radius of the sphere for cart_geo='geo'.

Details

This function contains three different methods for trajectory calculations on unstructured triangular meshes for two-dimensional velocity fields: the well-known (first-order) explicit Euler method, the second-order Petterssen scheme, and a piece-wise analytical approach, Pansy. The latter is yet only available for Cartesian coordinates, but will be refurbished soon. All methods operate with an adaptive time step, as this both ensures efficiency of the algorithm and reduces errors.

There are a lot of input parameters, but most of them can usually be ignored. For large numbers of elements within the grid, I recommend calculating i.neighs and tri.cont.ini once in advance for repeated use. As a rule of thumb, the default time step size should roughly equal (or be lower than) the ratio of a common element circumradius and a representative value for speed.

The time grid is a numerical vector here, not a 'time' object. Therefore you need to ensure manually that the input time grid tgr (and eventually the radius Rsphere) matches the units of the velocity vectors. If U is for instance given in meter per seconds, it is one possible solution to

• set Rsphere to a value in [m], for instance 6371e3

• set tgr to values in [s], for instance 0:10*86400 for daily velocities

• set dt.def to a value in [s], for instance 0.5*86400 for half a day.

This also holds for T.end and t.ini.

Value

A list containing

pos	a numerical array (dim: [N.ini, 3, N.time.out.max]) containing the trajectories ([,1:2,]) and the elements of the given position ([,3,]) for the output time grid of length N.time.out.max, which is set by maxiter, dt.def and qfrac.
time	a numerical array (dim: [N.ini, N.time.out.max]) containing the individual output time grids for each trajectory. Due to the adaptive time step sizes, these do not necessarily coincide.
flag	a vector of length N.ini which eventually contains flags for each trajectory under certain conditions (<will be revised soon>). If all entries are NA, it is very likely that the trajectories meet your expectations.

Author(s)

Simon Reifenberg

References

<work in progress>

Examples

#################################################################
############### Example for Cartesian coordinates ###############
#################################################################

# generate and load example data (exchange filepath 'fp' at will)
fp = "~/"
sl.tracer.generate.testcase(filepath = fp)
load(paste0(fp, "sl.tracer.example.data.RData"))        # contains U, x.ini, grid, i.neighs, tri.cont.ini

##### TRAJECTORY COMPUTATION ####
dt.def  = 1     # set default time step size, for playing around I recommend values between 0.1 and 4

# only 'grid', 'U', 'x.ini' are used from the example data to keep the function call comparably tidy here
# -> T.end, t.ini, tgr, maxiter, qfrac are set to default, tri.cont.ini and i.neights are calculated during runtime
trajPe <- sl.tracer.calculate.trajectories(U, x.ini, grid=grid, cart_geo="cart", dt.def = dt.def, method="Petterssen")
trajEu <- sl.tracer.calculate.trajectories(U, x.ini, grid=grid, cart_geo="cart", dt.def = dt.def, method="Euler")
trajPa <- sl.tracer.calculate.trajectories(U, x.ini, grid=grid, cart_geo="cart", dt.def = dt.def, method="Pansy")

##### PLOT RESULTS #####
#(the analytical solution would be on concentric circles, black dots mark initial positions)
par(pty="s", mfrow = c(1,3))
# Petterssen
pir <- sl.plot.init("lonlat", lonlat.lonrange = c(-7,7), lonlat.latrange = c(-7,7), do.init.device = F)
sl.plot.elem(pir, grid$lon, grid$lat, grid$elem, col.border = "grey", fill = F, lwd = .5)
for (itraj in 1:10){sl.plot.lines(pir, lon = trajPe$pos[itraj,1,], lat = trajPe$pos[itraj,2,], col = "purple", lwd  = 1.5)}
sl.plot.points(pir, lon = x.ini[1,], lat = x.ini[2,], col = "black", pch = 19, cex = .7) # init. positions
sl.plot.end(pir, do.close.device = F); mtext("Petterssen", col = "gray40", cex = 1.2)

# Euler
pir <- sl.plot.init("lonlat", lonlat.lonrange = c(-7,7), lonlat.latrange = c(-7,7), do.init.device = F)
sl.plot.elem(pir, grid$lon, grid$lat, grid$elem, col.border = "grey", fill = F, lwd = .5)
for (itraj in 1:10){sl.plot.lines(pir, lon = trajEu$pos[itraj,1,], lat = trajEu$pos[itraj,2,], col = "steelblue", lwd  = 1.5)}
sl.plot.points(pir, lon = x.ini[1,], lat = x.ini[2,], col = "black", pch = 19, cex = .7) # init. positions
sl.plot.end(pir, do.close.device = F); mtext("Euler", col = "gray40", cex = 1.2)

# Pansy
pir <- sl.plot.init("lonlat", lonlat.lonrange = c(-7,7), lonlat.latrange = c(-7,7), do.init.device = F)
sl.plot.elem(pir, grid$lon, grid$lat, grid$elem, col.border = "grey", fill = F, lwd = .5)
for (itraj in 1:10){sl.plot.lines(pir, lon = trajPa$pos[itraj,1,], lat = trajPa$pos[itraj,2,], col = "orange", lwd  = 1.5)}
sl.plot.points(pir, lon = x.ini[1,], lat = x.ini[2,], col = "black", pch = 19, cex = .7) # init. positions
mtext("Pansy", col = "gray40", cex = 1.2); sl.plot.end(pir, do.close.device = F)

FESOM/spheRlab documentation built on April 6, 2024, 6:52 p.m.