docs + parallel tracelines

DESCRIPTION

@@ -15,6 +15,7 @@ BugReports: https://github.com/cneyens/raem/issues
 Imports: 
     deSolve,
     graphics,
+    parallel,
     stats
 Suggests: 
     knitr,

R/aem.R

@@ -365,9 +365,9 @@ resfac <- function(element, aem) {
   return(resfac)
 }
 
-#' Add or remove an element to existing `aem` object
+#' Add or remove an element from an existing `aem` object
 #'
-#' [add_element()] adds a new element to the `aem` object.
+#' [add_element()] adds a new element to or from an `aem` object.
 #'
 #' @param aem `aem` object
 #' @param element analytic element of class `element`

R/tracelines.R

@@ -149,6 +149,7 @@ outside_vertical <- function(aem, x, y, z, ...) {
 #' @param R numeric, retardation coefficient passed to [velocity()]. Defaults to 1 (no retardation).
 #' @param tfunc function or list of functions with additional termination events for particles. See details. Defaults to `NULL`.
 #' @param tol numeric tolerance used to define when particles have crossed a line element. Defaults to 0.001 length units.
+#' @param ncores integer, number of cores to use when running in parallel. Defaults to 0 (no parallel computing). See details.
 #' @param ... ignored
 #'
 #' @details [deSolve::lsoda()] is used to numerically integrate the velocity vector.
@@ -171,6 +172,9 @@ outside_vertical <- function(aem, x, y, z, ...) {
 #'
 #' Backward particle tracking is performed by reversing the flow field (i.e. multiplying the velocities with `-1`).
 #'
+#' Traceline computation is embarrassingly parallel. When `ncores > 0`, the `parallel` package is used to set up the cluster with the requested nodes and
+#'     the tracelines are computed using [parallel::parLapplyLB()]. `ncores` should not exceed the number of available cores as returned by [parallel::detectCores()].
+#'
 #' @return [tracelines()] returns an object of class `tracelines` which is a list with length equal to the number of particles where each list element contains
 #'    a matrix with columns `time`, `x`, `y` and `z` specifying the registered time and coordinates of the particle as is it tracked through the flow field.
 #'
@@ -252,7 +256,15 @@ outside_vertical <- function(aem, x, y, z, ...) {
 #' # plot vertical cross-section of traceline 4 along increasing y-axis (from south to north)
 #' plot(paths[[4]][,c('y', 'z')], type = 'l')
 #'
-tracelines <- function(aem, x0, y0, z0, times, forward = TRUE, R = 1, tfunc = NULL, tol = 1e-3, ...) {
+#' @examplesIf parallel::detectCores() > 1
+#' # parallel computing
+#' m <- aem(k, top, base, n = n, uf, rf)
+#'
+#' x0 <- -200; y0 <- seq(-200, 200, 50)
+#' times <- seq(0, 25 * 365, 365 / 4)
+#' paths <- tracelines(m, x0 = x0, y0 = y0, z = top, times = times, ncores = 2)
+#'
+tracelines <- function(aem, x0, y0, z0, times, forward = TRUE, R = 1, tfunc = NULL, tol = 1e-3, ncores = 0, ...) {
 
   if(!is.null(tfunc) & !is.list(tfunc)) tfunc <- list(tfunc)
   direction <- ifelse(forward, 1, -1)
@@ -305,8 +317,35 @@ tracelines <- function(aem, x0, y0, z0, times, forward = TRUE, R = 1, tfunc = NU
                    rootfun = rootfun)
   }, SIMPLIFY = FALSE)
 
-  # get paths and clean
-  paths <- get_paths(x0, y0, z0)
+
+  # get paths
+  if(ncores > 0) { # parallel
+    if(ncores > parallel::detectCores()) stop('ncores > available cores', call. = FALSE)
+
+    crds <- data.frame(x = x0, y = y0, z = z0)
+    crds_list <- unname(split(crds, seq(nrow(crds))))
+
+    # make cluster
+    clust <- parallel::makeCluster(ncores)
+
+    try({
+      load_raem <- parallel::clusterCall(clust, function(...) library(raem))
+
+      parallel::clusterExport(clust,
+                  list("aem", "R", "vxvy", "rootfun", "get_paths",
+                       "times", "direction", "outside_vertical", "tfunc"),
+                  envir = environment())
+
+      paths <- parallel::parLapplyLB(clust, crds_list, function(i) get_paths(i[[1]], i[[2]], i[[3]])[[1]])
+    })
+
+    parallel::stopCluster(clust)
+
+  } else { # sequential
+    paths <- get_paths(x0, y0, z0)
+  }
+
+  # clean
   paths.m <- lapply(paths, matrix, ncol = 4, dimnames = list(NULL, c('time', 'x', 'y', 'z')))
   names(paths.m) <- NULL
   # if(length(paths) == 1) paths.m <- paths.m[[1]] # return matrix instead of list of matrices when only one particle is tracked

tests/testthat/test-tracelines.R

@@ -211,3 +211,32 @@ test_that('initial particles are not stuck', {
   )
 
 })
+
+test_that("parallel tracelines work", {
+
+  ncores <- 2 # CRAN only allows 2 cores
+  if(parallel::detectCores() < ncores) skip(paste('ncores < ', ncores))
+
+  k <- 10
+  top <- 10; base <- 0
+  n <- 0.2
+  R <- 1.5
+  i <- 0.001
+  alpha <- -45
+  TR <- k * (top - base)
+  hc <- 15
+
+  uf <- uniformflow(TR = TR, gradient = i, angle = alpha)
+  rf <- constant(-1000, 0, hc = hc)
+  m <- aem(k, top, base, n = n, uf, rf)
+
+  x0 <- seq(-200, 200, length = 10)
+  y0 <- 0
+  times <- seq(0, 365, 365 / 20)
+
+  paths <- tracelines(m, x0 = x0, y0 = y0, z = top, times = times, R = R)
+  pathsp <- tracelines(m, x0 = x0, y0 = y0, z = top, times = times, R = R, ncores = ncores)
+
+  expect_equal(paths, pathsp)
+
+})

vignettes/overview.Rmd

@@ -254,7 +254,7 @@ The use of `add_element()` to set-up the model, which can be used to add element
 
 A major benefit of the AEM approach is that no numerical grid of the domain is needed. This allows to zoom in on an area of interest without altering or rerunning the model, as illustrated by the plot of the inset area around well A.
 
-```{r example-model, fig.show="hold", out.width='50%'}
+```{r example-model, fig.show="hold"}
 
 # aquifer parameters ----
 k = 15          # hydraulic conductivity, m/d
@@ -336,6 +336,10 @@ for(i in m$elements) {
 
 # Output
 
+We'll use the model created above to show the different types of output that are available in `raem`.
+
 # Particle traces
 
+## Capture zone
+
 # Additional functionality