View source: R/sl.tracer.calculate.trajectories.R
sl.tracer.calculate.trajectories | R Documentation |
A suite of three different integration methods for kinematic trajectory calculations on unstructured triangular meshes.
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)
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: |
x.ini |
a numerical array (dim: |
t.ini |
a numerical scalar, the time for which the initial positions |
grid |
a list containing at least the entries |
i.neighs |
a list of length |
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 |
dt.def |
a numeric scalar defining the default time step size (in units of |
T.end |
a numeric scalar defining the maximal length of the time integration (in units of |
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 |
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
.
A list containing
pos |
a numerical array (dim: |
time |
a numerical array (dim: |
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 |
Simon Reifenberg
<work in progress>
#################################################################
############### 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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.