[TF24] Prototype of multirate stepper (DO NOT MERGE)

R/RcppExports.R

@@ -4089,6 +4089,18 @@ TF24_Environment__compute_rates <- function(obj_, resource_depletion) {
     invisible(.Call('_plant_TF24_Environment__compute_rates', PACKAGE = 'plant', obj_, resource_depletion))
 }
 
+TF24_Environment__r_analytic_partial_flow <- function(obj_, theta, dt) {
+    .Call('_plant_TF24_Environment__r_analytic_partial_flow', PACKAGE = 'plant', obj_, theta, dt)
+}
+
+TF24_Environment__r_residual_rhs <- function(obj_, theta, resource_depletion) {
+    .Call('_plant_TF24_Environment__r_residual_rhs', PACKAGE = 'plant', obj_, theta, resource_depletion)
+}
+
+TF24_Environment__r_drainage_touchdown_time <- function(obj_, theta, layer) {
+    .Call('_plant_TF24_Environment__r_drainage_touchdown_time', PACKAGE = 'plant', obj_, theta, layer)
+}
+
 TF24_Environment__time__get <- function(obj_) {
     .Call('_plant_TF24_Environment__time__get', PACKAGE = 'plant', obj_)
 }
@@ -4213,6 +4225,22 @@ make_node_schedule__Parameters___K93__K93_Env <- function(p) {
     .Call('_plant_make_node_schedule__Parameters___K93__K93_Env', PACKAGE = 'plant', p)
 }
 
+mri_fast_rate_calls_get <- function() {
+    .Call('_plant_mri_fast_rate_calls_get', PACKAGE = 'plant')
+}
+
+mri_fast_rate_calls_reset <- function() {
+    invisible(.Call('_plant_mri_fast_rate_calls_reset', PACKAGE = 'plant'))
+}
+
+patch_rhs_calls_get <- function() {
+    .Call('_plant_patch_rhs_calls_get', PACKAGE = 'plant')
+}
+
+patch_rhs_calls_reset <- function() {
+    invisible(.Call('_plant_patch_rhs_calls_reset', PACKAGE = 'plant'))
+}
+
 #' Generate a suitable set of default node introduction times,
 #' biased so that introductions are more closely packed at the
 #' beginning of time, become increasingly spread out.

R/RcppR6.R

@@ -1,6 +1,6 @@
 ## Generated by RcppR6: do not edit by hand
 ## Version: 0.2.4
-## Hash: e40f3d10382824a03e692ce606540515
+## Hash: 423758920dc13c34159566f729c5e6b5
 
 ##' @importFrom Rcpp evalCpp
 ##' @importFrom R6 R6Class
@@ -5627,6 +5627,15 @@ StochasticPatchRunner <- function(T, E) {
       },
       compute_rates = function(resource_depletion) {
         TF24_Environment__compute_rates(self, resource_depletion)
+      },
+      r_analytic_partial_flow = function(theta, dt) {
+        TF24_Environment__r_analytic_partial_flow(self, theta, dt)
+      },
+      r_residual_rhs = function(theta, resource_depletion) {
+        TF24_Environment__r_residual_rhs(self, theta, resource_depletion)
+      },
+      r_drainage_touchdown_time = function(theta, layer) {
+        TF24_Environment__r_drainage_touchdown_time(self, theta, layer)
       }),
     active=list(
       time = function(value) {

inst/RcppR6_classes.yml

@@ -295,6 +295,9 @@ Control:
   list:
     - function_integration_rule: size_t
     - shading_model: std::string
+    - ode_method: std::string
+    - n_collocation_nodes: size_t
+    - mri_use_split: bool
     - ppa_layer_optical_depth: double
     - ppa_layer_smoothing: double
     - offspring_production_tol: double
@@ -1359,6 +1362,15 @@ TF24_Environment:
     compute_rates:
       args: [resource_depletion: std::vector<double>]
       return_type: void
+    r_analytic_partial_flow:
+      args: [theta: std::vector<double>, dt: double]
+      return_type: std::vector<double>
+    r_residual_rhs:
+      args: [theta: std::vector<double>, resource_depletion: std::vector<double>]
+      return_type: std::vector<double>
+    r_drainage_touchdown_time:
+      args: [theta: double, layer: int]
+      return_type: double
   active:
     time: {type: double, access: field}
     soil_moist_sat: {type: double, access: field}

inst/include/plant/RcppR6_post.hpp

@@ -315,6 +315,9 @@ template <> inline SEXP wrap(const plant::Control& x) {
   Rcpp::List ret;
   ret["function_integration_rule"] = Rcpp::wrap(x.function_integration_rule);
   ret["shading_model"] = Rcpp::wrap(x.shading_model);
+  ret["ode_method"] = Rcpp::wrap(x.ode_method);
+  ret["n_collocation_nodes"] = Rcpp::wrap(x.n_collocation_nodes);
+  ret["mri_use_split"] = Rcpp::wrap(x.mri_use_split);
   ret["ppa_layer_optical_depth"] = Rcpp::wrap(x.ppa_layer_optical_depth);
   ret["ppa_layer_smoothing"] = Rcpp::wrap(x.ppa_layer_smoothing);
   ret["offspring_production_tol"] = Rcpp::wrap(x.offspring_production_tol);
@@ -355,6 +358,12 @@ template <> inline plant::Control as(SEXP x) {
   ret.function_integration_rule = Rcpp::as<size_t >(xl["function_integration_rule"]);
   // ret.shading_model = Rcpp::as<decltype(retshading_model) >(xl["shading_model"]);
   ret.shading_model = Rcpp::as<std::string >(xl["shading_model"]);
+  // ret.ode_method = Rcpp::as<decltype(retode_method) >(xl["ode_method"]);
+  ret.ode_method = Rcpp::as<std::string >(xl["ode_method"]);
+  // ret.n_collocation_nodes = Rcpp::as<decltype(retn_collocation_nodes) >(xl["n_collocation_nodes"]);
+  ret.n_collocation_nodes = Rcpp::as<size_t >(xl["n_collocation_nodes"]);
+  // ret.mri_use_split = Rcpp::as<decltype(retmri_use_split) >(xl["mri_use_split"]);
+  ret.mri_use_split = Rcpp::as<bool >(xl["mri_use_split"]);
   // ret.ppa_layer_optical_depth = Rcpp::as<decltype(retppa_layer_optical_depth) >(xl["ppa_layer_optical_depth"]);
   ret.ppa_layer_optical_depth = Rcpp::as<double >(xl["ppa_layer_optical_depth"]);
   // ret.ppa_layer_smoothing = Rcpp::as<decltype(retppa_layer_smoothing) >(xl["ppa_layer_smoothing"]);

inst/include/plant/control.h

@@ -33,6 +33,34 @@ struct Control {
   // so it never costs a string comparison on the hot path.
   std::string shading_model;
 
+  // Number of collocation nodes for the multirate (method="mri") fast sub-cycle.
+  // The soil sub-cycle needs the per-layer root uptake at each micro-step; that
+  // uptake is a density-weighted integral of per-cohort consumption over the
+  // size distribution. 0 (the default) evaluates it over all N cohorts (exact);
+  // m > 0 quadratures it at m of the N frozen cohorts instead, so the sub-cycle
+  // costs m physiology solves, not N. Converges ~O(m^-2), but the accuracy at a
+  // given m depends strongly on the distribution: young stands reach the sub-1%
+  // range by m≈20, mature (skewed) stands need m close to N for the same
+  // accuracy under the current even-index node placement (a smarter,
+  // importance-weighted placement is the open item -- plant#53 §6/§4.4). Only
+  // consulted on the method="mri" path; ignored otherwise.
+  size_t n_collocation_nodes;
+
+  // Multirate (method="mri") fast-block inner stepper. false (default) sub-cycles
+  // the soil block with the adaptive black-box RK; true uses the exact-flow split
+  // (R1 analytic drainage recession + ROS34PW2 on the gentle remainder), which
+  // removes the drainage stiffness so the sub-cycle takes far fewer micro steps
+  // (Lever 1). Only consulted on the method="mri" path.
+  bool mri_use_split;
+
+  // ODE integration method for the SCM resident solver. One of "rkck" (the
+  // default adaptive Cash-Karp explicit RK), "rodas" (the stiff Rosenbrock
+  // stepper), or "mri" (the multirate MRI-GARK stepper: a fixed macro grid that
+  // sub-cycles the fast soil column, for TF24). Empty is treated as "rkck", so
+  // default behaviour is unchanged. Selected once when the SCM builds its
+  // Solver; every other integration path is untouched.
+  std::string ode_method;
+
   // PPA only: thickness of one discrete canopy layer, in optical-depth units
   // (tau = sum of k * leaf-area-index above a height). The stepped light
   // profile floors tau to integer multiples of this value. The default 0.5

inst/include/plant/individual.h

@@ -144,6 +144,23 @@ template <typename T, typename E> class Individual {
     return rate(HEIGHT_INDEX);
   }
 
+  // Per-layer resource consumption of an individual grown to `height` under the
+  // given (frozen-light, soil-θ) environment. The factored coupling the
+  // multirate collocation consumes (plant#53 item 4): the stand's per-layer
+  // uptake is a density-weighted integral of this over the size distribution, so
+  // it can be quadratured at m << N heights instead of every cohort. Mirrors
+  // growth_rate_given_height: set the height, recompute, read the result.
+  std::vector<double> consumption_given_height(double height,
+                                               const environment_type& environment) {
+    set_state(HEIGHT_INDEX, height);
+    compute_rates(environment);
+    std::vector<double> ret(environment.ode_size());
+    for (size_t i = 0; i < ret.size(); ++i) {
+      ret[i] = consumption_rate(i);
+    }
+    return ret;
+  }
+
   double resource_compensation_point() {
     environment_type env = environment_type();
 

inst/include/plant/models/tf24_environment.h

@@ -281,11 +281,7 @@ class TF24_Environment : public Environment {
 
     double water_input;
     double rainfall = extrinsic_drivers.evaluate("rainfall", time);
-    const double soil_moist_sat_0 =
-      soil_parameter_value(soil_moist_sat_layers, soil_moist_sat, 0);
-    double infiltration = rainfall * std::max(
-      0.0,
-      1 - a_infil * std::pow(vars.state(0) / soil_moist_sat_0, b_infil));
+    double infiltration = infiltration_rate(vars.state(0));
     double total_resource_depletion = 0;
 
 
@@ -335,6 +331,103 @@ class TF24_Environment : public Environment {
 
   }
 
+  // ------------------------------------------------------------------
+  // R1 operator split of the soil water balance
+  // ------------------------------------------------------------------
+  // The single-layer gravitational drainage is a power-law loss with a
+  // closed-form recession, so it is integrated exactly; infiltration, the
+  // inter-layer cascade inflow and root uptake are the gentle remainder the
+  // caller steps. Composed (flow, remainder, flow) they reproduce compute_rates.
+  // This removes the wet-end drainage stiffness (and its positivity clamp) from
+  // any soil integration, single-rate included, and is the exact-flow half an
+  // operator-splitting stepper needs. These operate on the soil moisture only;
+  // the cumulative-flux diagnostics stay with compute_rates.
+
+  // Infiltration into the top layer: rainfall reduced by saturation-excess
+  // runoff. One source of truth for compute_rates and residual_rhs.
+  double infiltration_rate(double theta0) const {
+    const double rainfall = extrinsic_drivers.evaluate("rainfall", time);
+    const double sat0 = soil_parameter_value(soil_moist_sat_layers, soil_moist_sat, 0);
+    return rainfall * std::max(0.0, 1 - a_infil * std::pow(theta0 / sat0, b_infil));
+  }
+
+  // Advance the per-layer drainage theta' = -K(theta)/dz exactly over dt, in
+  // place. K(theta) = K_sat (theta/theta_sat)^p, p = 2 n_psi + 3, so
+  // theta' = -c theta^p with c = K_sat/(dz theta_sat^p) and the recession
+  //   theta(dt) = [theta^{1-p} + (p-1) c dt]^{-1/(p-1)}
+  // is positivity-preserving. Above saturation K caps at K_sat (a linear decline)
+  // until theta reaches theta_sat, matching soil_K_from_soil_theta's clamp.
+  void analytic_partial_flow(std::vector<double>& theta, double dt) const {
+    for (size_t i = 0; i < soil_number_of_depths; i++) {
+      const double k_sat = soil_parameter_value(K_sat_layers, K_sat, i);
+      const double sat = soil_parameter_value(soil_moist_sat_layers, soil_moist_sat, i);
+      const double p = 2 * soil_parameter_value(n_psi_layers, n_psi, i) + 3;
+      double th = theta[i], rem = dt;
+      // Drainage stops at the residual floor (issue #485): a layer is not drained
+      // below theta_r. This is the exact-flow form of the monolithic positivity
+      // guard, and matches drainage_touchdown_time.
+      if (th <= soil_moist_residual) continue;
+      if (th > sat) {                                // capped drainage: linear decline
+        const double t_to_sat = (th - sat) * dz[i] / k_sat;
+        if (rem <= t_to_sat) { theta[i] = th - k_sat * rem / dz[i]; continue; }
+        th = sat; rem -= t_to_sat;
+      }
+      const double c = k_sat / (dz[i] * std::pow(sat, p));
+      const double drained =
+        std::pow(std::pow(th, 1 - p) + (p - 1) * c * rem, -1.0 / (p - 1));
+      theta[i] = std::max(drained, soil_moist_residual);
+    }
+  }
+
+  // The non-drainage remainder: infiltration into layer 0, the drainage cascade
+  // inflow K(theta_{i-1}) into deeper layers, and root uptake. drate[i] is the
+  // rate the caller steps between exact-flow half-steps. The uptake floor (issue
+  // #485) keeps a layer at the residual moisture from being dried further by
+  // uptake; drainage positivity is already guaranteed by analytic_partial_flow.
+  void residual_rhs(const std::vector<double>& theta,
+                    const std::vector<double>& resource_depletion,
+                    std::vector<double>& drate) const {
+    const double infiltration = infiltration_rate(theta[0]);
+    for (size_t i = 0; i < soil_number_of_depths; i++) {
+      const double water_in =
+        (i == 0) ? infiltration : soil_K_from_soil_theta(theta[i - 1], i - 1);
+      double rate = (water_in - resource_depletion[i]) / dz[i];
+      if (theta[i] <= soil_moist_residual && !(rate > 0.0)) rate = 0.0;
+      drate[i] = rate;
+    }
+  }
+
+  // Time for the drainage recession to bring a layer from theta down to the
+  // residual floor -- a closed-form contact event for a splitting micro-stepper,
+  // rather than a dense-output root-find. Infinity if already at/below the floor.
+  double drainage_touchdown_time(double theta, size_t layer) const {
+    if (theta <= soil_moist_residual) return std::numeric_limits<double>::infinity();
+    const double sat = soil_parameter_value(soil_moist_sat_layers, soil_moist_sat, layer);
+    const double k_sat = soil_parameter_value(K_sat_layers, K_sat, layer);
+    const double p = 2 * soil_parameter_value(n_psi_layers, n_psi, layer) + 3;
+    double t = 0.0, th = theta;
+    if (th > sat) { t += (th - sat) * dz[layer] / k_sat; th = sat; }
+    const double c = k_sat / (dz[layer] * std::pow(sat, p));
+    t += (std::pow(soil_moist_residual, 1 - p) - std::pow(th, 1 - p)) / ((p - 1) * c);
+    return t;
+  }
+
+  // Value-in/value-out wrappers so the split can be exercised and checked against
+  // the monolithic solve from R (which cannot observe an in-place update).
+  std::vector<double> r_analytic_partial_flow(std::vector<double> theta, double dt) const {
+    analytic_partial_flow(theta, dt);
+    return theta;
+  }
+  std::vector<double> r_residual_rhs(std::vector<double> theta,
+                                     std::vector<double> resource_depletion) const {
+    std::vector<double> drate(soil_number_of_depths);
+    residual_rhs(theta, resource_depletion, drate);
+    return drate;
+  }
+  double r_drainage_touchdown_time(double theta, int layer) const {
+    return drainage_touchdown_time(theta, static_cast<size_t>(layer));
+  }
+
   // calculate K from K_sat based on theta
   double soil_K_from_soil_theta(double theta, size_t layer) const {
     //Eq. 5 Zeng and Decker (2009), ref Clapp and Hornberger (1978)