Implementation

implementation

This page overviews the numerical methods used in plant and serves as the index for the Implementation section, linking out to detailed pages for each method.

Characteristic method

System of ODEs, IVPs.

See more detail in the characteristic method.

ODE stepping

See more detail in ODE stepping and control.

Adaptive node spacing

See more details in the node spacing algorithm.

TK splines

Gauss-Kronrod quadrature

To numerically integrate functions, we use the Gauss-Kronrod quadrature method. This is implemented in the QK and QAG classes.

  • General background on the method is available at: https://en.wikipedia.org/wiki/Gaussian_quadrature
  • The QK class is used to numerically integrate functions using Gauss-Kronrod quadrature.
    • The number of points used when numerically integrating the function is defined by the control variable function_integration_rule. The integration has several “rules”, including QK15, QK21, QK31, QK41, QK51, defined in qk_rules.cpp. These determine the numbers of points in the integration. The default number of points is 21.
    • The QK code was ported from GSL by Rich FitzJohn.
  • The QAG class uses the QK class, for adaptive integration using Gauss-Kronrod quadrature. In the adaptive method, the difference between a Gauss quadrature rule and its Kronrod extension is used as an estimate of the approximation error, from which the need for more points is determined.
    • The QAG routine is adapted from the “QAG” algorithm in QUADPACK.
    • The QAG class calls class QK, which does the actual integration.

Numerical derivatives

node introduction

Trapezium integration

trapezium(xx, yy).

Pre-compute step

Approximating the light environment with a spline

Creating the light environment

The patch calls compute_environment,

  • which calls environment.compute_environment(f, height_max()), where f = compute_competition
    • compute_competition passes through species
      • in species, compute_competition integrates density over nodes
        • in node, compute_competition = density * individual.compute_competition(height_);
          • individual passes straight through to strategy
            • strategy calculates fraction of leaf above = k_I * area_leaf(height) * Q(z, height);
  • environment.compute_environment creates a spline of light availability, as a function of height
    • to do this, it uses the ResourceSpline class.
      • ResourceSpline class uses an adaptive routine to sample across the height range to approximate the function
        • the function is ultimately approximated with a spline, using the tk class
    • creating the spline is achieved through calls to compute_competition, as described above

For more on light environments, see canopy models.

Using the light environment
  • In compute_rates, environment gets passed from Patch -> Species -> Node -> Individual -> Strategy
    • FF16_strategy net_mass_production_dt calls assimilation
      • assimilation uses the object quadrature::QK function_integrator to integrate the function assimilation_leaf(environment.get_environment_at_height(z)) * q(z, height) over z, where z is depth in the plant. The function is integrated using Gauss-Kronrod quadrature. The number of points used when numerically integrating the function is defined by the argument function_integration_rule (see above), the default number of points is 21.
        • Previously, we used the QAG (instead of QK) class, to enable adaptive integration of the function. The QAG class uses QK, but includes an extra layer to refine the number of points. However, this is slower and isn’t noticeably changing the results. The adaptive routine may take twice as long.