Warn if w_max too small

NEWS.md

@@ -1,5 +1,7 @@
 # development version
 
+- Adds a warning to `project()` and `projectToSteady()` / `steady()` if the abundance of any species at its maximum size `w_max` is not negligible, indicating that the chosen `w_max` is too small and individuals are being cut off by the upper boundary condition.
+
 - New experimental `steadyNewton()` finds a steady state by solving the
   steady-state equation directly with a Newton-type root finder (using the
   `nleqslv` package) instead of running the dynamics to convergence. Unlike

R/helpers.R

@@ -302,3 +302,52 @@ get_steady_state_n <- function(params, g, mu, D, N0,
 
     return(zero_above_support(n, w_top))
 }
+
+#' Warn if w_max is too small for any species
+#'
+#' Checks whether the total flux (advective + diffusive) of any species at its
+#' maximum size `w_max` is non-negligible compared to the maximum total flux
+#' across all sizes.
+#' 
+#' It also issues a warning if species do not even reach maturity size.
+#'
+#' @param params A MizerParams object.
+#' @param n An array (species x size) with the number density.
+#' @param rates A list of rates as returned by `mizerRates()`.
+#' @param threshold The threshold ratio above which a warning is issued.
+#' @noRd
+warn_if_w_max_violation <- function(params, n, rates, threshold = 1e-5) {
+    w_top <- support_top_idx(params)
+    no_sp <- nrow(params@species_params)
+    flux <- getFlux(params, n = n)
+    violating_species <- character()
+    starving_species <- character()
+    for (i in seq_len(no_sp)) {
+        w_mat_idx <- sum(params@w <= params@species_params[i, "w_mat"])
+        w_max_idx <- sum(params@w <= params@species_params[i, "w_max"])
+        # increase by 1 is needed in second-order scheme and it is also
+        # needed in first-order scheme to not warn on existing models
+        w_max_idx <- min(w_max_idx + 1, length(params@w))
+        if (flux[i, w_mat_idx] > 0) {
+            ratio <- flux[i, w_top[i]] / flux[i, w_mat_idx]
+            if (ratio > threshold) {
+                violating_species <- c(violating_species, params@species_params$species[i])
+            }
+        } else {
+            starving_species <- c(starving_species, params@species_params$species[i])
+        }
+    }
+    if (length(violating_species) > 0) {
+        warning("For the following species, the abundance is not negligible at their maximum size w_max:\n",
+                paste(violating_species, collapse = ", "),
+                ".\nThis means that individuals are being cut off by the upper boundary condition. ",
+                "You should choose a larger w_max for these species.",
+                call. = FALSE)
+    }
+    if (length(starving_species) > 0) {
+        warning("The following species never reach their maturity size:\n",
+                paste(violating_species, collapse = ", "),
+                call. = FALSE)
+    }
+}
+

R/project.R

@@ -391,6 +391,8 @@ project.MizerParams <- function(object, effort,
         sim@n_other[i, ] <- unserialize(serialize(n_list$n_other, NULL))
     }
 
+    warn_if_w_max_violation(params, n = n_list$n, rates = n_list$rates)
+
     return(sim)
 }
 

R/steady.R

@@ -218,6 +218,8 @@ projectToSteady.MizerParams <- function(params,
     params@initial_n_pp[] <- current$n_pp
     params@initial_n_other[] <- current$n_other
 
+    warn_if_w_max_violation(params, n = current$n, rates = current$rates)
+
     if (return_sim) {
         sim@params <- params
         sel <- 1:i

tests/testthat/test-project.R

@@ -498,3 +498,15 @@ test_that("Predation diffusion changes simulation trajectory", {
     sim_base <- project(NS_params_small, t_max = 1)
     expect_false(identical(sim_d@n, sim_base@n))
 })
+
+test_that("warn_if_w_max_violation warns when w_max is too small", {
+    # Check that it doesn't warn on default NS_params_small
+    expect_warning(project(NS_params_small, t_max = 1), NA)
+
+    # Check that it warns when w_max of Herring is set too small
+    params <- NS_params_small
+    params@species_params$w_max[2] <- 5
+    params <- validParams(params)
+    expect_warning(project(params, t_max = 1), "abundance is not negligible at their maximum size w_max")
+})
+

tests/testthat/test-steadyNewton.R

@@ -100,15 +100,17 @@ test_that("steadyNewton works with the second-order (van Leer) scheme", {
     sow$flux <- "van_leer"
     sow$bin_average <- TRUE
     second_order_w(p) <- sow
-    ps <- steady(p, t_max = 100, progress_bar = FALSE)
+    expect_warning(ps <- steady(p, t_max = 100, progress_bar = FALSE),
+                   "abundance is not negligible")
 
     pn <- steadyNewton(ps)
     expect_s4_class(pn, "MizerParams")
     # The returned state must still be a fixed point of the (second-order)
     # dynamics. The van Leer limiter is only Lipschitz, so the Newton residual
     # cannot be driven to machine precision, but the projected drift is the
     # honest test and stays small.
-    sim <- project(pn, t_max = 1, dt = 0.25, t_save = 1)
+    expect_warning(sim <- project(pn, t_max = 1, dt = 0.25, t_save = 1),
+                   "abundance is not negligible")
     n0 <- pn@initial_n
     n1 <- finalN(sim)
     support <- n0 > 0