diff --git a/R/RcppExports.R b/R/RcppExports.R index 8420932d..53cef3d1 100644 --- a/R/RcppExports.R +++ b/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. diff --git a/R/RcppR6.R b/R/RcppR6.R index e043ffc3..38b7e6d3 100644 --- a/R/RcppR6.R +++ b/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) { diff --git a/inst/RcppR6_classes.yml b/inst/RcppR6_classes.yml index e11cb572..0cb949f5 100644 --- a/inst/RcppR6_classes.yml +++ b/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] return_type: void + r_analytic_partial_flow: + args: [theta: std::vector, dt: double] + return_type: std::vector + r_residual_rhs: + args: [theta: std::vector, resource_depletion: std::vector] + return_type: std::vector + 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} diff --git a/inst/include/plant/RcppR6_post.hpp b/inst/include/plant/RcppR6_post.hpp index 91fdecaf..ce43f251 100644 --- a/inst/include/plant/RcppR6_post.hpp +++ b/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(xl["function_integration_rule"]); // ret.shading_model = Rcpp::as(xl["shading_model"]); ret.shading_model = Rcpp::as(xl["shading_model"]); + // ret.ode_method = Rcpp::as(xl["ode_method"]); + ret.ode_method = Rcpp::as(xl["ode_method"]); + // ret.n_collocation_nodes = Rcpp::as(xl["n_collocation_nodes"]); + ret.n_collocation_nodes = Rcpp::as(xl["n_collocation_nodes"]); + // ret.mri_use_split = Rcpp::as(xl["mri_use_split"]); + ret.mri_use_split = Rcpp::as(xl["mri_use_split"]); // ret.ppa_layer_optical_depth = Rcpp::as(xl["ppa_layer_optical_depth"]); ret.ppa_layer_optical_depth = Rcpp::as(xl["ppa_layer_optical_depth"]); // ret.ppa_layer_smoothing = Rcpp::as(xl["ppa_layer_smoothing"]); diff --git a/inst/include/plant/control.h b/inst/include/plant/control.h index 95d5fdfe..757c0061 100644 --- a/inst/include/plant/control.h +++ b/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 diff --git a/inst/include/plant/individual.h b/inst/include/plant/individual.h index 94939ee0..c6ee56e6 100644 --- a/inst/include/plant/individual.h +++ b/inst/include/plant/individual.h @@ -144,6 +144,23 @@ template 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 consumption_given_height(double height, + const environment_type& environment) { + set_state(HEIGHT_INDEX, height); + compute_rates(environment); + std::vector 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(); diff --git a/inst/include/plant/models/tf24_environment.h b/inst/include/plant/models/tf24_environment.h index 059b793c..6f7522ec 100644 --- a/inst/include/plant/models/tf24_environment.h +++ b/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& 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& theta, + const std::vector& resource_depletion, + std::vector& 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::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 r_analytic_partial_flow(std::vector theta, double dt) const { + analytic_partial_flow(theta, dt); + return theta; + } + std::vector r_residual_rhs(std::vector theta, + std::vector resource_depletion) const { + std::vector 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(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) diff --git a/inst/include/plant/patch.h b/inst/include/plant/patch.h index c83a7062..9e7b7185 100644 --- a/inst/include/plant/patch.h +++ b/inst/include/plant/patch.h @@ -6,6 +6,11 @@ #include #include #include +// Patch declares an MRI fast/slow partition (below), so it is a multirate +// System: pull in the multirate stepper's out-of-line definition here, so every +// TU that builds a Solver has MriStep::step defined and not just +// declared. Harmless for method != "mri" (the stepper is never invoked). +#include #include @@ -16,6 +21,21 @@ using namespace Rcpp; namespace plant { +// Diagnostic counter (defined in src/mri_diag.cpp): number of fast-block +// coupling evaluations across a run, the multirate cost driver. Incremented in +// fast_rates; read/reset from R via mri_fast_rate_calls_get/reset. +extern long mri_fast_rate_calls; +extern long patch_rhs_calls; + +// Does the environment expose the R1 operator split (exact drainage recession + +// residual)? Only TF24 does; the multirate split inner (Lever 1) is gated on it +// so a patch over an environment without it (e.g. FF16) still compiles. +template +struct env_has_split : std::false_type {}; +template +struct env_has_split().analytic_partial_flow( + std::declval&>(), 0.0))>> : std::true_type {}; + template class Patch { public: @@ -96,6 +116,118 @@ class Patch { // Retrieve auxillary variables and save into the ode solver odelia::ode::iterator ode_aux(odelia::ode::iterator it) const; + // * Multirate (MRI) partition interface + // Additive hooks that let odelia's method="mri" stepper treat the patch as a + // fast/slow system. They are only touched by MriStep; every other integration + // path is untouched, so production runs are bit-identical when method != "mri". + // + // The ODE state is already laid out [cohorts | environment] (see ode_state), + // which is exactly MRI's [slow | fast] layout: the cohorts are the slow block + // and the soil column is the fast block, contiguous at the tail. No coupling + // aggregate is used -- the fast block reads its slow context from the frozen + // light field captured once per leg by freeze_slow -- so coupling_size is 0. + size_t slow_size() const { return ode_size() - environment.ode_size(); } + size_t fast_size() const { return environment.ode_size(); } + size_t coupling_size() const { return 0; } + void aggregate(const std::vector&, std::vector&) const {} + + // Freeze the slow (cohort) context for one MRI leg: set the cohort states and + // rebuild the light field they cast, so the ensuing fast sub-cycle varies only + // the soil column against a fixed canopy. + void freeze_slow(const std::vector& x) { + odelia::ode::set_ode_state(species.begin(), species.end(), x.begin()); + compute_environment(true); + environment_ptr = &environment; + // Freeze the collocation nodes against this leg's distribution (see + // fast_rates); a no-op when collocation is off. + if (control.n_collocation_nodes > 0) { + for (auto& s : species) { + s.set_collocation_nodes(control.n_collocation_nodes); + } + } + } + + // Per-layer root uptake at the current (soil-θ) environment: the density- + // weighted cohort consumption / area. Full-N re-solves every cohort's + // physiology; collocation quadratures it over m frozen cohorts. Shared by the + // adaptive and split fast-block hooks. + std::vector fast_block_uptake() { + if (control.n_collocation_nodes > 0) { + std::vector depletion(environment.ode_size(), 0.0); + for (auto& s : species) { + s.add_collocated_consumption(environment, depletion); + } + for (auto& d : depletion) { + d /= area; + } + return depletion; + } + compute_species_rates(); + return assemble_resource_depletion(); + } + + // Fast tendency: the soil (environment) rates at soil state u, with the canopy + // frozen by the preceding freeze_slow. Uptake is re-evaluated at this soil + // moisture (it depends on θ), then the environment rates are read out. g is + // unused (coupling_size == 0). + void fast_rates(const std::vector& u, const std::vector& /*g*/, + std::vector& du) { + ++mri_fast_rate_calls; + environment.set_ode_state(u.begin()); + environment_ptr = &environment; + environment.compute_rates(fast_block_uptake()); + environment.ode_rates(du.begin()); + } + + // Multirate split inner (Lever 1): opt in via control when the environment + // provides the R1 exact-flow hooks. Off (adaptive inner) otherwise. + bool mri_split() const { + return env_has_split::value && control.mri_use_split; + } + + // Exact stiff-drainage flow on the soil layers of the fast block (aux slots -- + // the cumulative-flux diagnostics -- are advanced by the residual, not here). + void analytic_flow(std::vector& u, double dt) { + if constexpr (env_has_split::value) { + const size_t ns = static_cast(environment.get_soil_number_of_depths()); + std::vector theta(u.begin(), u.begin() + ns); + environment.analytic_partial_flow(theta, dt); + std::copy(theta.begin(), theta.end(), u.begin()); + } + } + + // Gentle remainder the split inner steps: infiltration + inter-layer cascade − + // uptake on the soil layers (drainage is handled exactly by analytic_flow). g + // is unused (coupling_size == 0); uptake is re-evaluated at this soil state. + // The trailing aux (cumulative-flux) slots are left at zero rate here -- they + // are diagnostics decoupled from the dynamics. + void residual_rhs(const std::vector& u, const std::vector& /*g*/, + std::vector& du) { + ++mri_fast_rate_calls; + std::fill(du.begin(), du.end(), 0.0); + if constexpr (env_has_split::value) { + environment.set_ode_state(u.begin()); + environment_ptr = &environment; + const std::vector depletion = fast_block_uptake(); + const size_t ns = static_cast(environment.get_soil_number_of_depths()); + std::vector theta(u.begin(), u.begin() + ns), drate(ns); + environment.residual_rhs(theta, depletion, drate); + std::copy(drate.begin(), drate.end(), du.begin()); + } + } + + // Slow tendency: the cohort rates at cohort state x and soil state u. The + // canopy has moved, so the light field is rebuilt before the rates are read. + void slow_rates(const std::vector& x, const std::vector& u, + std::vector& dx) { + odelia::ode::set_ode_state(species.begin(), species.end(), x.begin()); + environment.set_ode_state(u.begin()); + compute_environment(true); + environment_ptr = &environment; + compute_rates(); + odelia::ode::ode_rates(species.begin(), species.end(), dx.begin()); + } + // Returns state in structure format as opposed to single // vector as given by ode_state Rcpp::List r_get_state() const; @@ -164,6 +296,9 @@ class Patch { int idx = 0; // used to access environment cache for mutant runs void compute_environment(bool rescale); void compute_rates(); + // Reusable pieces of compute_rates (shared with the multirate split hooks). + void compute_species_rates(); + std::vector assemble_resource_depletion() const; // Seed the patch from parameters.initial_state (nodes + birth bookkeeping) // when present; called from reset(). Sets environment.time = initial_time. @@ -184,9 +319,6 @@ class Patch { environment_type environment; std::vector species; - //TODO(#476): Move into environment? - std::vector resource_depletion; - environment_type* environment_ptr; Control control; @@ -258,9 +390,6 @@ void Patch::reset() { s.resize_consumption_rates(environment.ode_size()); } - // resize to species count - resource_depletion.reserve(environment.ode_size()); - // compute ephemeral effects like light_availability environment.clear(); @@ -573,37 +702,43 @@ void Patch::compute_environment(bool rescale) { template void Patch::compute_rates() { - - // Computes rates of change for the patch, including all the component species + ++patch_rhs_calls; + // Computes rates of change for the patch, including all the component species. // While the patch has an `environment`, the rates here are calculated from - // the env_ptr, which is a pointer to an environment object - // -- for the resident the pointer points to the internal environment object - // -- for a mutant, the pointer points to a cached environment object - double time_ = environment_ptr->time; + // environment_ptr: + // -- for the resident it points to the internal environment object + // -- for a mutant it points to a cached environment object + compute_species_rates(); + environment_ptr->compute_rates(assemble_resource_depletion()); +} +// The per-species rate evaluation (growth, mortality, fecundity, and the +// per-cohort resource consumption read back below). Split out of compute_rates +// so the multirate fast-block hooks can reuse it without the monolithic soil +// rate assembly. +template +void Patch::compute_species_rates() { + double time_ = environment_ptr->time; double pr_patch_survival = survival_weighting->pr_survival(time_); for (size_t i = 0; i < size(); ++i) { double birth_rate = species[i].extrinsic_drivers().evaluate("birth_rate", time_); - - // Pass the environment that pointer is tracking into compute rates. species[i].compute_rates(*environment_ptr, pr_patch_survival, birth_rate); } +} - resource_depletion.reserve(environment_ptr->ode_size()); - for(size_t i = 0; i < environment_ptr->ode_size(); i++) { - double resource_consumed = std::accumulate(species.begin(), species.end(), 0.0, [i](double r, const species_type& s) { - return r + s.consumption_rate(i); // accumulates r from zero - }); - - resource_depletion.push_back(resource_consumed/area); +// Per-resource depletion = density-weighted cohort consumption / area, one entry +// per environment ODE slot (soil layers used by the environment rates; the +// trailing aux slots carry no consumption). Requires compute_species_rates first. +template +std::vector Patch::assemble_resource_depletion() const { + std::vector depletion; + depletion.reserve(environment_ptr->ode_size()); + for (size_t i = 0; i < environment_ptr->ode_size(); i++) { + double resource_consumed = std::accumulate(species.begin(), species.end(), 0.0, + [i](double r, const species_type& s) { return r + s.consumption_rate(i); }); + depletion.push_back(resource_consumed / area); } - - - environment_ptr->compute_rates(resource_depletion); - - //todo do we need to clear this every step? - resource_depletion.clear(); - + return depletion; } // TODO(#478): We should only be recomputing the light environment for the diff --git a/inst/include/plant/scm.h b/inst/include/plant/scm.h index b8c6a416..8d343f5b 100644 --- a/inst/include/plant/scm.h +++ b/inst/include/plant/scm.h @@ -15,6 +15,20 @@ using namespace Rcpp; namespace plant { +// Map the Control's ode_method string onto the odelia Solver's Method enum. +// Empty is treated as "rkck" so the default is unchanged; an unknown value is +// an error rather than a silent fall-through to the default. +inline odelia::ode::Method scm_ode_method(const Control& control) { + const std::string& m = control.ode_method; + if (m.empty() || m == "rkck") return odelia::ode::Method::rkck; + if (m == "rodas") return odelia::ode::Method::rodas; + if (m == "imex") return odelia::ode::Method::imex; + if (m == "mri") return odelia::ode::Method::mri; + util::stop("Unknown control$ode_method '" + m + + "'; use 'rkck', 'rodas' or 'mri'."); + return odelia::ode::Method::rkck; // unreachable; silences the return warning +} + // SCM: the "Solver for Characteristics Method" driver. // // Owns a Patch (the population being integrated), a NodeSchedule (when each @@ -139,7 +153,7 @@ template SCM::SCM(parameters_type p, environment_type e, Control c) : parameters(p), control(c), patch(parameters, e, c), node_schedule(make_node_schedule(parameters)), - solver(patch, make_ode_control(c)) { + solver(patch, make_ode_control(c), scm_ode_method(c)) { parameters.validate(); diff --git a/inst/include/plant/species.h b/inst/include/plant/species.h index ff83d6fb..7cd1d8ec 100644 --- a/inst/include/plant/species.h +++ b/inst/include/plant/species.h @@ -69,6 +69,21 @@ class Species : public SpeciesBase, T, E, Node> { double consumption_rate(int i) const; std::vector consumption_rate_by_node_rev(int i) const; + // * Collocation of the per-layer uptake (plant#53 item 4) + // consumption_rate(i) is a density-weighted trapezium of per-cohort uptake + // over the size distribution. When the canopy (and hence the distribution) is + // frozen -- e.g. across a multirate leg -- that integral can be evaluated at + // m << N heights instead of every node. set_collocation_nodes captures m + // log-spaced heights spanning the current distribution and their + // density*trapezium weights; add_collocated_consumption then re-evaluates only + // those m cohorts under the current (soil-θ) environment and adds their + // density-weighted uptake into `depletion`. Convergence is O(m^-2) (measured: + // <0.5% by m≈15-20). A no-op when fewer than two nodes exist (matching + // consumption_rate, which needs a distribution to integrate). + void set_collocation_nodes(size_t m); + void add_collocated_consumption(const environment_type& environment, + std::vector& depletion); + odelia::ode::iterator ode_aux(odelia::ode::iterator it) const; Rcpp::NumericMatrix r_get_state() const; @@ -125,6 +140,14 @@ class Species : public SpeciesBase, T, E, Node> { using base_type::control; node_type new_node; + // Collocation scratch (see set_collocation_nodes / add_collocated_consumption). + // colloc_probe_ is a single reusable individual grown to each node height in + // turn; colloc_weights_[n] folds the density at colloc_heights_[n] into the + // trapezium weight, so the aggregate is one weighted sum over m re-solves. + individual_type colloc_probe_; + std::vector colloc_heights_; + std::vector colloc_weights_; + typedef typename std::vector::iterator nodes_iterator; typedef typename std::vector::const_iterator nodes_const_iterator; }; @@ -132,7 +155,8 @@ class Species : public SpeciesBase, T, E, Node> { template Species::Species(strategy_type s) : base_type(s), - new_node(this->strategy) { + new_node(this->strategy), + colloc_probe_(this->strategy) { } template @@ -300,6 +324,54 @@ std::vector Species::consumption_rate_by_node_rev(int i) const { return ret; } +template +void Species::set_collocation_nodes(size_t m) { + colloc_heights_.clear(); + colloc_weights_.clear(); + const size_t n = size(); + if (n < 2 || m < 2) { + return; // no distribution to integrate (consumption_rate returns 0 too) + } + + // Subsample m of the N nodes (evenly by index, endpoints included) and + // quadrature over their *actual* heights and densities -- the same trapezium + // consumption_rate uses, on a coarser node set. Using real (height, density) + // pairs avoids interpolating the irregular size distribution onto synthetic + // heights, which is the dominant collocation error. If m >= N we just take + // every node (exact). Nodes are stored descending; walk them into ascending + // height order for the trapezium. + const size_t mm = std::min(m, n); + std::vector h(mm), density(mm); + for (size_t j = 0; j < mm; ++j) { + // even index in [0, n-1], ascending height (node n-1 is the shortest). + const size_t idx = (mm == 1) ? 0 : (j * (n - 1)) / (mm - 1); + const node_type& node = nodes[n - 1 - idx]; + h[j] = node.height(); + density[j] = node.get_density(); + } + + colloc_heights_ = h; + colloc_weights_.resize(mm); + for (size_t j = 0; j < mm; ++j) { + const double w = (j == 0) ? 0.5 * (h[1] - h[0]) + : (j == mm - 1) ? 0.5 * (h[mm - 1] - h[mm - 2]) + : 0.5 * (h[j + 1] - h[j - 1]); + colloc_weights_[j] = w * density[j]; + } +} + +template +void Species::add_collocated_consumption(const environment_type& environment, + std::vector& depletion) { + for (size_t j = 0; j < colloc_heights_.size(); ++j) { + const std::vector c = + colloc_probe_.consumption_given_height(colloc_heights_[j], environment); + for (size_t i = 0; i < depletion.size(); ++i) { + depletion[i] += c[i] * colloc_weights_[j]; + } + } +} + // bit clunky... template size_t Species::aux_size() const { diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 586a3526..51897403 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -11503,6 +11503,45 @@ BEGIN_RCPP return R_NilValue; END_RCPP } +// TF24_Environment__r_analytic_partial_flow +std::vector TF24_Environment__r_analytic_partial_flow(plant::RcppR6::RcppR6 obj_, std::vector theta, double dt); +RcppExport SEXP _plant_TF24_Environment__r_analytic_partial_flow(SEXP obj_SEXP, SEXP thetaSEXP, SEXP dtSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< plant::RcppR6::RcppR6 >::type obj_(obj_SEXP); + Rcpp::traits::input_parameter< std::vector >::type theta(thetaSEXP); + Rcpp::traits::input_parameter< double >::type dt(dtSEXP); + rcpp_result_gen = Rcpp::wrap(TF24_Environment__r_analytic_partial_flow(obj_, theta, dt)); + return rcpp_result_gen; +END_RCPP +} +// TF24_Environment__r_residual_rhs +std::vector TF24_Environment__r_residual_rhs(plant::RcppR6::RcppR6 obj_, std::vector theta, std::vector resource_depletion); +RcppExport SEXP _plant_TF24_Environment__r_residual_rhs(SEXP obj_SEXP, SEXP thetaSEXP, SEXP resource_depletionSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< plant::RcppR6::RcppR6 >::type obj_(obj_SEXP); + Rcpp::traits::input_parameter< std::vector >::type theta(thetaSEXP); + Rcpp::traits::input_parameter< std::vector >::type resource_depletion(resource_depletionSEXP); + rcpp_result_gen = Rcpp::wrap(TF24_Environment__r_residual_rhs(obj_, theta, resource_depletion)); + return rcpp_result_gen; +END_RCPP +} +// TF24_Environment__r_drainage_touchdown_time +double TF24_Environment__r_drainage_touchdown_time(plant::RcppR6::RcppR6 obj_, double theta, int layer); +RcppExport SEXP _plant_TF24_Environment__r_drainage_touchdown_time(SEXP obj_SEXP, SEXP thetaSEXP, SEXP layerSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< plant::RcppR6::RcppR6 >::type obj_(obj_SEXP); + Rcpp::traits::input_parameter< double >::type theta(thetaSEXP); + Rcpp::traits::input_parameter< int >::type layer(layerSEXP); + rcpp_result_gen = Rcpp::wrap(TF24_Environment__r_drainage_touchdown_time(obj_, theta, layer)); + return rcpp_result_gen; +END_RCPP +} // TF24_Environment__time__get double TF24_Environment__time__get(plant::RcppR6::RcppR6 obj_); RcppExport SEXP _plant_TF24_Environment__time__get(SEXP obj_SEXP) { @@ -11850,6 +11889,44 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// mri_fast_rate_calls_get +double mri_fast_rate_calls_get(); +RcppExport SEXP _plant_mri_fast_rate_calls_get() { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + rcpp_result_gen = Rcpp::wrap(mri_fast_rate_calls_get()); + return rcpp_result_gen; +END_RCPP +} +// mri_fast_rate_calls_reset +void mri_fast_rate_calls_reset(); +RcppExport SEXP _plant_mri_fast_rate_calls_reset() { +BEGIN_RCPP + Rcpp::RNGScope rcpp_rngScope_gen; + mri_fast_rate_calls_reset(); + return R_NilValue; +END_RCPP +} +// patch_rhs_calls_get +double patch_rhs_calls_get(); +RcppExport SEXP _plant_patch_rhs_calls_get() { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + rcpp_result_gen = Rcpp::wrap(patch_rhs_calls_get()); + return rcpp_result_gen; +END_RCPP +} +// patch_rhs_calls_reset +void patch_rhs_calls_reset(); +RcppExport SEXP _plant_patch_rhs_calls_reset() { +BEGIN_RCPP + Rcpp::RNGScope rcpp_rngScope_gen; + patch_rhs_calls_reset(); + return R_NilValue; +END_RCPP +} // node_schedule_times_default std::vector node_schedule_times_default(double max_time); RcppExport SEXP _plant_node_schedule_times_default(SEXP max_timeSEXP) { @@ -13043,6 +13120,9 @@ static const R_CallMethodDef CallEntries[] = { {"_plant_TF24_Environment__get_atm_o2_kpa", (DL_FUNC) &_plant_TF24_Environment__get_atm_o2_kpa, 1}, {"_plant_TF24_Environment__get_atm_kpa", (DL_FUNC) &_plant_TF24_Environment__get_atm_kpa, 1}, {"_plant_TF24_Environment__compute_rates", (DL_FUNC) &_plant_TF24_Environment__compute_rates, 2}, + {"_plant_TF24_Environment__r_analytic_partial_flow", (DL_FUNC) &_plant_TF24_Environment__r_analytic_partial_flow, 3}, + {"_plant_TF24_Environment__r_residual_rhs", (DL_FUNC) &_plant_TF24_Environment__r_residual_rhs, 3}, + {"_plant_TF24_Environment__r_drainage_touchdown_time", (DL_FUNC) &_plant_TF24_Environment__r_drainage_touchdown_time, 3}, {"_plant_TF24_Environment__time__get", (DL_FUNC) &_plant_TF24_Environment__time__get, 1}, {"_plant_TF24_Environment__time__set", (DL_FUNC) &_plant_TF24_Environment__time__set, 2}, {"_plant_TF24_Environment__soil_moist_sat__get", (DL_FUNC) &_plant_TF24_Environment__soil_moist_sat__get, 1}, @@ -13074,6 +13154,10 @@ static const R_CallMethodDef CallEntries[] = { {"_plant_TF24f_oderunner_individual_internals", (DL_FUNC) &_plant_TF24f_oderunner_individual_internals, 1}, {"_plant_node_schedule_default__Parameters___K93__K93_Env", (DL_FUNC) &_plant_node_schedule_default__Parameters___K93__K93_Env, 1}, {"_plant_make_node_schedule__Parameters___K93__K93_Env", (DL_FUNC) &_plant_make_node_schedule__Parameters___K93__K93_Env, 1}, + {"_plant_mri_fast_rate_calls_get", (DL_FUNC) &_plant_mri_fast_rate_calls_get, 0}, + {"_plant_mri_fast_rate_calls_reset", (DL_FUNC) &_plant_mri_fast_rate_calls_reset, 0}, + {"_plant_patch_rhs_calls_get", (DL_FUNC) &_plant_patch_rhs_calls_get, 0}, + {"_plant_patch_rhs_calls_reset", (DL_FUNC) &_plant_patch_rhs_calls_reset, 0}, {"_plant_node_schedule_times_default", (DL_FUNC) &_plant_node_schedule_times_default, 1}, {"_plant_FF16_strategy_expand_allometry", (DL_FUNC) &_plant_FF16_strategy_expand_allometry, 4}, {"_plant_TF24_strategy_expand_allometry", (DL_FUNC) &_plant_TF24_strategy_expand_allometry, 4}, diff --git a/src/RcppR6.cpp b/src/RcppR6.cpp index d7d69903..91412b46 100644 --- a/src/RcppR6.cpp +++ b/src/RcppR6.cpp @@ -4628,6 +4628,18 @@ void TF24_Environment__compute_rates(plant::RcppR6::RcppR6compute_rates(resource_depletion); } // [[Rcpp::export]] +std::vector TF24_Environment__r_analytic_partial_flow(plant::RcppR6::RcppR6 obj_, std::vector theta, double dt) { + return obj_->r_analytic_partial_flow(theta, dt); +} +// [[Rcpp::export]] +std::vector TF24_Environment__r_residual_rhs(plant::RcppR6::RcppR6 obj_, std::vector theta, std::vector resource_depletion) { + return obj_->r_residual_rhs(theta, resource_depletion); +} +// [[Rcpp::export]] +double TF24_Environment__r_drainage_touchdown_time(plant::RcppR6::RcppR6 obj_, double theta, int layer) { + return obj_->r_drainage_touchdown_time(theta, layer); +} +// [[Rcpp::export]] double TF24_Environment__time__get(plant::RcppR6::RcppR6 obj_) { return obj_->time; } diff --git a/src/control.cpp b/src/control.cpp index 697df5e8..9d80ce08 100644 --- a/src/control.cpp +++ b/src/control.cpp @@ -19,6 +19,18 @@ Control::Control() { // unchanged for both. shading_model = ""; + // ODE method for the SCM resident solver (see Control header). Empty = "rkck", + // so default behaviour is unchanged. + ode_method = ""; + + // Multirate fast sub-cycle collocation nodes (see Control header). 0 = exact + // full-N uptake; only used on the method="mri" path. + n_collocation_nodes = 0; + + // Multirate fast-block inner stepper (see Control header). false = adaptive + // black-box RK; true = exact-flow split. Only used on the method="mri" path. + mri_use_split = false; + // PPA canopy layer thickness in optical-depth units (see Control header). ppa_layer_optical_depth = 0.5; // PPA layer-boundary smoothing fraction (see Control header). diff --git a/src/leaf_model.cpp b/src/leaf_model.cpp index 11070b94..baa597bb 100644 --- a/src/leaf_model.cpp +++ b/src/leaf_model.cpp @@ -671,6 +671,13 @@ void Leaf::set_shutdown_state(double root_collar) { root_collar_psi_ = root_collar; opt_psi_stem_ = psi_crit; profit_ = -R_d_ - hydraulic_cost_TF(psi_crit); + // The stem is held at psi_crit -- transpiration is not possible -- so the + // plant draws no water. Without this the reused leaf keeps whatever + // soil_consumption_/E_up_ the previous (responsive) cohort left, feeding a + // stale, soil-moisture-independent phantom uptake into the soil balance and + // zeroing the reverse-mode gradient across the drought regime (R-C, #55). + soil_consumption_.assign(soil_number_of_depths_, 0.0); + E_up_ = 0.0; } // Shared setup + feasibility handling for the root-collar solve. Extracted diff --git a/src/mri_diag.cpp b/src/mri_diag.cpp new file mode 100644 index 00000000..66c9f9f3 --- /dev/null +++ b/src/mri_diag.cpp @@ -0,0 +1,38 @@ +// Diagnostic counter for the multirate (method="mri") fast sub-cycle. +// +// Each Patch::fast_rates call evaluates the soil (fast-block) coupling once, +// which costs O(M) cohort physiology solves -- the fast-block cost driver. This +// counter accumulates those calls across a run so we can read the per-macro-step +// fast-evaluation count (Lever-1 headroom: how many micro evaluations the stiff +// drainage currently forces) and compare schemes. Not part of the model; used +// only by the multirate cost measurements. +#include + +namespace plant { +long mri_fast_rate_calls = 0; +// Full patch RHS evaluations (Patch::compute_rates), i.e. one evaluation of the +// whole coupled derivative -- the cost/step proxy shared by every global stepper +// (rkck, rodas). Lets us compare accepted-step economics across methods and +// diagnose whether the coupled system is accuracy- or stability-limited. +long patch_rhs_calls = 0; +} + +// [[Rcpp::export]] +double mri_fast_rate_calls_get() { + return static_cast(plant::mri_fast_rate_calls); +} + +// [[Rcpp::export]] +void mri_fast_rate_calls_reset() { + plant::mri_fast_rate_calls = 0; +} + +// [[Rcpp::export]] +double patch_rhs_calls_get() { + return static_cast(plant::patch_rhs_calls); +} + +// [[Rcpp::export]] +void patch_rhs_calls_reset() { + plant::patch_rhs_calls = 0; +} diff --git a/tests/testthat/test-control.R b/tests/testthat/test-control.R index a1e43369..fae4b225 100644 --- a/tests/testthat/test-control.R +++ b/tests/testthat/test-control.R @@ -17,6 +17,9 @@ test_that("Defaults", { function_integration_rule = 21, # size_t so not int shading_model = "", # empty = each strategy's own default + ode_method = "", # empty = "rkck" (the SCM resident solver default) + n_collocation_nodes = 0, # size_t; 0 = exact full-N multirate uptake + mri_use_split = FALSE, # multirate split inner off by default ppa_layer_optical_depth = 0.5, ppa_layer_smoothing = 0.3, diff --git a/tests/testthat/test-tf24-shutdown.R b/tests/testthat/test-tf24-shutdown.R new file mode 100644 index 00000000..93c334dd --- /dev/null +++ b/tests/testthat/test-tf24-shutdown.R @@ -0,0 +1,58 @@ +# R-C / #55: a shut-down TF24 leaf (stem held at psi_crit, transpiration impossible) +# draws no water. Before the fix, set_shutdown_state left soil_consumption_/E_up_ +# stale, so a reused leaf fed a soil-moisture-independent phantom uptake into the +# soil balance and zeroed the reverse-mode gradient across the drought regime. + +mkleaf <- function() Leaf(vcmax_25 = 100, jmax_25 = 100 * 167, c = 2.04, b = 3, psi_crit = 5, + root_c = 2.65, root_b = 1.29, root_psi_crit = 1.29 * (log(1 / 0.05))^(1 / 2.65), beta2 = 1, + a = 0.3, curv_fact_elec_trans = 0.7, curv_fact_colim = 0.99, GSS_tol_abs = 1e-8, + vulnerability_curve_ncontrol = 100, ci_abs_tol = 1e-8, ci_niter = 1000, g1_TF24 = 46.32995, + beta_R_H = 3.4e3, beta_R_V = 9.4e4) +setp <- function(L, psi) L$set_physiology(area_leaf = 1, mass_root_prop = rep(1 / length(psi), length(psi)), + rho = 608, a_bio = 0.0245, PPFD = 1800, psi_soil = psi, soil_depth = seq(0.3, by = 0.3, length.out = length(psi)), + leaf_specific_conductance_max = 1e-4, atm_vpd = 1, ca = 40, sapwood_volume_per_leaf_area = 1e-4, + leaf_temp = 25, atm_o2_kpa = 21, atm_kpa = 101.3) + +test_that("a shut-down TF24 leaf draws no water (#55)", { + L <- mkleaf() + setp(L, rep(6, 5)) # every layer drier than psi_crit (= 5) -> shutdown + L$find_root_collar_psi() + expect_true(is.finite(L$profit_)) + expect_identical(L$E_up_, 0) + expect_identical(L$soil_consumption_, rep(0, 5)) +}) + +test_that("a reused leaf does not carry stale uptake into shutdown (#55)", { + L <- mkleaf() + setp(L, rep(1, 5)); L$find_root_collar_psi() # responsive solve leaves nonzero soil_consumption_ + expect_gt(sum(L$soil_consumption_), 0) + setp(L, rep(6, 5)); L$find_root_collar_psi() # shutdown on the SAME leaf must reset to 0 + expect_identical(L$E_up_, 0) + expect_identical(L$soil_consumption_, rep(0, 5)) +}) + +test_that("TF24 stand uptake vanishes in deep drought without a phantom floor (#55)", { + SAT <- 0.428; KSAT <- 163.0411; DZmm <- 300; P <- 2 * 6.57 + 3 + Kf <- function(th) KSAT * (pmax(th, 1e-9) / SAT)^P + p0 <- scm_base_parameters("TF24"); p0$max_patch_lifetime <- 30 + p1 <- add_strategies(p0, trait_matrix(0.0825, "lma")) + env <- Environment("TF24"); env$extrinsic_drivers_set_variable("rainfall", 0:2, c(0, 0, 0)) + h <- c(14, 10, 7, 5, 3.5, 2.2, 1.3, 0.7); d <- c(0.015, 0.03, 0.06, 0.12, 0.25, 0.6, 1.5, 4.0) + st <- make_initial_state(p1, heights = h, densities = d, env = env, ctrl = Control()) + patch <- SCM("TF24", "TF24_Env")(set_initial_state(p1, st), env, Control())$patch + patch$compute_environment() + ne <- patch$environment$ode_size; ns <- patch$ode_size; si <- (ns - ne + 1):(ns - ne + 5) + # aggregate per-layer uptake backed out of the real patch derivs at uniform theta (zero inflow) + agg_uptake <- function(theta) { + y <- patch$ode_state; y[si] <- pmin(pmax(rep(theta, 5), 1e-6), SAT - 1e-6) + r <- patch$derivs(y, 0)[si]; K <- Kf(theta) + up <- numeric(5); up[1] <- -K - DZmm / 1000 * r[1] # dz in m for this backing-out is immaterial to sign/zero + for (i in 2:5) up[i] <- -DZmm / 1000 * r[i] + up + } + expect_gt(sum(abs(agg_uptake(0.20))), 1e-4) # responsive when wet + expect_lt(max(abs(agg_uptake(0.04))), 1e-9) # ~0 in deep drought -- no phantom floor + # the drought-response gradient is live in the responsive regime (was dead before the fix) + slope <- (sum(agg_uptake(0.17)) - sum(agg_uptake(0.15))) / 0.02 + expect_gt(abs(slope), 1e-3) +}) diff --git a/tests/testthat/test-tf24-soil-split.R b/tests/testthat/test-tf24-soil-split.R new file mode 100644 index 00000000..90d5bf50 --- /dev/null +++ b/tests/testthat/test-tf24-soil-split.R @@ -0,0 +1,135 @@ +# R1 operator split of the TF24 soil water balance. The stiff per-layer drainage +# is integrated exactly (analytic_partial_flow); infiltration, the inter-layer +# cascade and uptake are the gentle residual (residual_rhs). Together they must +# reproduce the monolithic soil solve. + +# monolithic soil rate, re-derived independently in R from the TF24 physics, as an +# uncontaminated reference: theta' = (water_in - drainage - uptake) / dz. +soil_defaults <- list(sat = 0.428, Ksat = 163.0411, npsi = 6.57, depth = 1.5, + a_infil = 1, b_infil = 8, theta_res = 1e-2, n = 5) + +soil_K <- function(th, p) p$Ksat * (min(max(th, 0), p$sat) / p$sat)^(2 * p$npsi + 3) + +mono_rate <- function(theta, rd, rain, p) { + dz <- p$depth / p$n + infil <- rain * max(0, 1 - p$a_infil * (theta[1] / p$sat)^p$b_infil) + wout <- vapply(theta, soil_K, numeric(1), p = p) + rate <- numeric(length(theta)) + for (i in seq_along(theta)) { + win <- if (i == 1) infil else wout[i - 1] + r <- (win - wout[i] - rd[i]) / dz + if (theta[i] <= p$theta_res && !(r > 0)) r <- 0 + rate[i] <- r + } + rate +} + +rk4 <- function(theta, rate_fn, dt, nsub) { + h <- dt / nsub + for (s in seq_len(nsub)) { + k1 <- rate_fn(theta); k2 <- rate_fn(theta + 0.5 * h * k1) + k3 <- rate_fn(theta + 0.5 * h * k2); k4 <- rate_fn(theta + h * k3) + theta <- theta + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4) + } + theta +} + +# The exact statement that the split IS the monolithic RHS: the drainage-flow rate +# (-K/dz) plus the residual equals the monolithic rate, everywhere in the interior. +testthat::test_that("the split decomposes the monolithic soil rate exactly", { + p <- soil_defaults + env <- Environment("TF24") + env$set_soil_number_of_depths(p$n) + dz <- p$depth / p$n + cases <- list( + list(theta = c(0.42, 0.38, 0.30, 0.22, 0.15), rd = c(3, 2, 1.5, 1, 0.5), rain = 0), + list(theta = c(0.40, 0.35, 0.30, 0.25, 0.20), rd = rep(1, p$n), rain = 20), + list(theta = c(0.25, 0.24, 0.23, 0.22, 0.21), rd = c(2, 2, 2, 2, 2), rain = 5) + ) + for (cs in cases) { + env$extrinsic_drivers_set_constant("rainfall", cs$rain) + residual <- env$r_residual_rhs(cs$theta, cs$rd) + flow_rate <- -vapply(cs$theta, soil_K, numeric(1), p = p) / dz + mono <- mono_rate(cs$theta, cs$rd, cs$rain, p) + expect_equal(residual + flow_rate, mono, tolerance = 1e-10) + } +}) + +testthat::test_that("the exact drainage flow matches a fine integration", { + p <- soil_defaults + env <- Environment("TF24") + env$set_soil_number_of_depths(p$n) + dz <- p$depth / p$n + for (th0 in c(0.42, 0.30, 0.15)) { + for (dt in c(0.5, 2.0)) { + exact <- env$r_analytic_partial_flow(rep(th0, p$n), dt) + fine <- rk4(rep(th0, p$n), function(th) -vapply(th, soil_K, numeric(1), p = p) / dz, + dt, nsub = 20000) + expect_lt(max(abs(exact - fine)), 1e-7) + expect_true(all(exact > 0)) # positivity by construction + } + } +}) + +testthat::test_that("the touchdown time brings a layer to the residual floor", { + p <- soil_defaults + env <- Environment("TF24") + env$set_soil_number_of_depths(p$n) + for (th0 in c(0.42, 0.25)) { + t_touch <- env$r_drainage_touchdown_time(th0, 0) + expect_true(is.finite(t_touch) && t_touch > 0) + at_floor <- env$r_analytic_partial_flow(rep(th0, p$n), t_touch) + expect_equal(at_floor[1], p$theta_res, tolerance = 1e-6) + } + expect_true(is.infinite(env$r_drainage_touchdown_time(soil_defaults$theta_res, 0))) +}) + +# Strang integration converges to the monolithic solve. Where the monolithic +# explicit solve is well resolved (drought, mild forcing, at the floor) the split +# matches it closely; in the monsoon regime the explicit reference is itself +# stability-limited (the reason R1 exists), so there we only require the split to +# converge as the step shrinks and to stay bounded and positive. +split_run <- function(env, theta0, rd, T_end, nstep) { + theta <- theta0 + h <- T_end / nstep + for (step in seq_len(nstep)) { + theta <- env$r_analytic_partial_flow(theta, h / 2) + theta <- rk4(theta, function(th) env$r_residual_rhs(th, rd), h, nsub = 4) + theta <- env$r_analytic_partial_flow(theta, h / 2) + } + theta +} + +testthat::test_that("Strang split matches the monolithic solve when it is resolved", { + p <- soil_defaults + env <- Environment("TF24") + env$set_soil_number_of_depths(p$n) + T_end <- 4 + scenarios <- list( + list(rain = 0, rd = rep(2.0, p$n), theta0 = rep(0.40, p$n)), # drought + list(rain = 1, rd = c(3, 2.5, 2, 1.5, 1), theta0 = rep(0.30, p$n)), # mild + list(rain = 5, rd = rep(6.0, p$n), theta0 = rep(0.08, p$n)) # near the floor + ) + for (s in scenarios) { + env$extrinsic_drivers_set_constant("rainfall", s$rain) + mono <- rk4(s$theta0, function(th) mono_rate(th, s$rd, s$rain, p), T_end, nsub = 16000) + e_coarse <- max(abs(split_run(env, s$theta0, s$rd, T_end, 100) - mono)) + e_fine <- max(abs(split_run(env, s$theta0, s$rd, T_end, 400) - mono)) + expect_lt(e_fine, e_coarse) # Strang splitting error -> 0 with h + expect_lt(e_fine, 1e-2) # and is small at a modest step + } +}) + +testthat::test_that("Strang split converges and stays bounded in the stiff monsoon", { + p <- soil_defaults + env <- Environment("TF24") + env$set_soil_number_of_depths(p$n) + env$extrinsic_drivers_set_constant("rainfall", 30) + T_end <- 4; rd <- rep(1.0, p$n); theta0 <- rep(0.20, p$n) + mono <- rk4(theta0, function(th) mono_rate(th, rd, 30, p), T_end, nsub = 16000) + e1 <- max(abs(split_run(env, theta0, rd, T_end, 400) - mono)) + e2 <- max(abs(split_run(env, theta0, rd, T_end, 1600) - mono)) + expect_lt(e2, e1) # converging toward the reference + fine <- split_run(env, theta0, rd, T_end, 1600) + expect_true(all(is.finite(fine) & fine > 0 & fine < p$sat + 0.05)) +})