Tutorial 3: 3D geometry

In this tutorial, we will more advanced functionality relating to the 3D structure of our plant, including the simulation of different soil compartments and ray tracing. As our package relies on VirtualPlantLab.jl for most structural modelling functionality, it will play a key role in this notebook.

Setup

using Pkg; Pkg.activate("../..")

using PlutoUI; TableOfContents()

using PlantModules

using ModelingToolkit, OrdinaryDiffEq, Plots

import Plots: plot, plot!

using VirtualPlantLab, CairoMakie

import VirtualPlantLab: Mesh

import VirtualPlantLab.PlantGeomPrimitives: Vec

using SkyDomes, PlantBiophysics, 
    PlantBiophysics.PlantMeteo, PlantBiophysics.PlantSimEngine

using DataInterpolations

Context

In the previous two tutorials, we simulated water dynamics for a plant given a certain, static structure. This tutorial, we will also consider the structural growth of a plant. As PlantModules.jl has no inherent structural modelling capabilities, we will instead how to integrate with the graph rewriting functionality provided by VirtualPlantLab.jl to achieve this. Specifically, we will couple our package's functionality to the rewriting process by making the rewrite rules depend on the functional status of the tree.

In addition, this tutorial considers two more geometry-related functionalities typically expected of FSPMs:

We simulate carbon dynamics in the leaves based on a ray tracer, a carbon assimilation model and a stomatal conductance model, all of which are available in VirtualPlantLab.jl.
We discretize the soil into multiple compartments to simulate water transport more realistically.

Structural definition

As we will be using the ray tracing functionality from VirtualPlantLab.jl, we recommend users to read the corresponding tutorial first, as well as the prior tutorials in the same series. In a nutshell, the extra steps that are required are:

Define information about the materials of the structural modules that will interact with light;
Specify this material in the geometry functions;
Create one or more light sources;
Add the soil surface to take into account reflected light;
Create and run the ray tracer.

For didactical purposes, we will base this tutorial directly on the series of tutorials mentioned above. Note that in VirtualPlantLab.jl's ray tracing tutorial, light interception is used to simulate carbon dynamics based on a discrete sink-source implentation. We will deviate here and instead use PlantModules.jl's functionality to model carbon dynamics coupled to water relations based on differential equations.

The tree structure used in the ray tracing tutorial has a lot of additional information for the discrete sink-source simulation of carbon dynamics. Considering we don't require this information, we will use the simpler tree structure from the first tutorial in the series.

Shoots

Structural modules

The tree consists of the classic set of internodes, nodes and leaves, as well as buds that grow into branches, corresponding nodes left behind by the buds, and meristems that grow into additional phytomers.

Some things of note for this tutorial:

All plant parts that are plotted need to have their dimensions defined;
All plant parts that intercept light in the ray tracing additionally need to have a material defined (see the VPL docs);
Plant parts that are part of a rewriting rule depending on the plant's functional status, should have those functional variables defined and be defined to be mutable.

Base.@kwdef struct Meristem <: VirtualPlantLab.Node end

Base.@kwdef mutable struct Bud <: VirtualPlantLab.Node
    D::Vector{Float64} = [0.5]
    W::Float64 = PlantModules.volume(PlantModules.Sphere(), [0.5])
end

Bud

Base.@kwdef struct Node <: VirtualPlantLab.Node end

Base.@kwdef struct BudNode <: VirtualPlantLab.Node end

Base.@kwdef mutable struct Internode <: VirtualPlantLab.Node
    D::Vector{Float64} = [0.5, 3.0] # newly grown internodes start small
    mat::Lambertian{1} = Lambertian(τ = 0.05, ρ = 0.1) # internodes require a material as they interact with light
end

Internode

Base.@kwdef mutable struct Leaf <: VirtualPlantLab.Node
    D::Vector{Float64} = [5, 3, 0.05]
    mat::Lambertian{1} = Lambertian(τ = 0.05, ρ = 0.1) # leaves also interact with light
    PAR_samples::Vector{Float64} = Float64[]
    PAR_func::Function = (x...) -> error("A leaf node exists with no defined PAR function.")
end

Leaf

Geometry functions

Define a feed! method for every structural module that needs a defined geometry, i.e. that needs to be plotted or interacts with light.

function VirtualPlantLab.feed!(turtle::Turtle, i::Internode, vars)
    # Rotate turtle around the head to implement elliptical phyllotaxis
    rh!(turtle, vars.phyllotaxis)
    # remember to divide by 100 to convert cm to m
    HollowCylinder!(
        turtle, length = i.D[2] / 100, height = i.D[1] / 100,
        width = i.D[1] / 100, move = true, colors = RGB(0.5,0.4,0.0), materials = i.mat
    )
    return nothing
end

function VirtualPlantLab.feed!(turtle::Turtle, l::Leaf, vars)
    ra!(turtle, -vars.leaf_angle)
    Rectangle!(turtle, length = l.D[1] / 100, width = l.D[2] / 100, move = false,
             colors = RGB(0.2, 0.6, 0.2), materials = l.mat)
    ra!(turtle, vars.leaf_angle)
    
    return nothing
end

function VirtualPlantLab.feed!(turtle::Turtle, b::BudNode, vars)
    ra!(turtle, -vars.branch_angle)
end

Structural growth rules

The tree grows by two rewriting steps:

Tree meristems grow into phytomers, occuring every rewrite step,
Buds grow branches, occuring with a probability proportional to the amount of phytomers to the apical meristem,

and one elongation step:

Every internode elongates for a set fraction of its current length.

Base.@kwdef struct treeparams # parameters used in rewrite rules
    growth::Float64 = 0.1
    W_burst::Float64 = 1.5 # water content at which bud is guaranteed to burst
    phyllotaxis::Float64 = 140.0
    leaf_angle::Float64 = 30.0
    branch_angle::Float64 = 45.0
end

treeparams

meristem_rule = Rule(
    Meristem,
    rhs = mer -> Node() + (Bud(), Leaf()) + Internode() + Meristem()
)

Rule replacing nodes of type Meristem without context capturing.

function prob_break(bud)
    prob = data(bud).W / graph_data(bud).W_burst
    return rand() < prob
end

prob_break (generic function with 1 method)

branch_rule = Rule(
    Bud,
    lhs = prob_break,
    rhs = bud -> BudNode() + Internode() + Meristem()
)

Rule replacing nodes of type Bud without context capturing.

shoot_graph = Graph(
    axiom = Internode(D = [0.5, 10.0]) + Node() + (Bud(), Leaf()) + Internode() + Meristem(),
    rules = (meristem_rule, branch_rule), data = treeparams()
);

render(Mesh(shoot_graph))

Roots

We model the roots as consisting of a single structural module. At every rewriting step, it either elongates or splits in two with a set probability. As we will have to connect the root segments to the correct soil compartment later on, it's important to track the positions of the root segments. We include positional information in our model by assigning a growing direction to all root segments. When new root segments are produced, their location follows this growing direction with a certain random perturbation.

Base.@kwdef struct rootparams
    split_prob::Float64 = 0.2
    Δdirection_σ::Float64 = 0.2
    Δdirection_σ_branch::Float64 = 0.5
end;

struct Root <: VirtualPlantLab.Node
    D::Vector{Float64}
    coords::Vector{Float64}
    direction::Vector{Float64}
end

function move(root, Δdirection_σ)
    new_D = [0.5, 3.0] # newly grown roots start small
    new_direction = data(root).direction + Δdirection_σ * randn(3) |>
        x -> x / sqrt(sum(x.^2)) # normalize to unit length
    new_coords = data(root).coords + new_D[2] * new_direction
    
    moved_root = Root(
        new_D,
        new_coords,
        new_direction
    )
    
    return moved_root
end

move (generic function with 1 method)

root_rule = Rule(
    Root,
    lhs = root -> !has_children(root),
    rhs = root -> (
        rand() < graph_data(root).split_prob ? 
        Root(data(root).D, data(root).coords, data(root).direction) + 
            (
                move(root, graph_data(root).Δdirection_σ_branch),
                move(root, graph_data(root).Δdirection_σ_branch)
            ) : 
        Root(data(root).D, data(root).coords, data(root).direction) + 
            move(root, graph_data(root).Δdirection_σ)
    )
)

Rule replacing nodes of type Root without context capturing.

root_graph = Graph(
    axiom = sum([Root([0.5, 10.0], [0.0, 0.0, i*-10.0], [0.0, 0.0, -1.0]) for i in 1:11]),
    rules = (root_rule), data = rootparams()
)

Dynamic graph with 11 nodes of types Main.var"workspace#25".Root and 1 rewriting rules.
Dynamic graph variables stored in struct of type rootparams

Note

Purely to demonstrate connecting with soil voxels in the next part, we start off with a long root.

plotstructure(root_graph)

Environment

As mentioned before, we will consider a soil consisting of multiple compartments. More specifically, we will discretize the soil into voxels, or cubes, of a given size. To connect each voxel to its non-diagonal neighbours, we can simply define the graph as a 3-dimensional array of soil nodes, which will get translated into our desired graph by PlantModules.jl. However, as Julia arrays are no typical graph format, they store no information about the id, and therefore it is required to manually define the id inside of our structural modules. Additionally, we will include the \(x\), \(y\) and \(z\) coordinates of the voxel's centre.

Base.@kwdef struct Soil <: VirtualPlantLab.Node 
    id::Integer
    x
    y
    z
end

Soil

vs = 100.0; # voxel size

soil_graph = [
    Soil(id = x + 3*y + 9*z; x, y, z) 
    for x in vs*(-1:1), y in vs*(-1:1), z in vs*(-0.5:-1:-2.5)
];

plotstructure(soil_graph)

The air could be similarly divided into compartments. For simplicity, however, we will again model it as a single node.

struct Air <: VirtualPlantLab.Node end

Ray tracing

We create the light sources corresponding to a realistic sky, using the functionality from SkyDomes.jl (from the from VirtualPlantLab.jl ecosystem). The ray tracing tutorial we base ourselves on defines an average sky for an entire day, which is possible because they model carbon dynamics on the day scale. We deviate here and create a sky corresponding with a single point in time, as we simulate plant growth continuously throughout the day.

waveband_conversion(Itype = :diffuse, waveband = :PAR, mode = :power)

0.6

function create_sky(day_fraction; mesh, lat = 52.0*π/180.0, DOY = 182)
    # Compute solar irradiance
    sky_light = clear_sky(lat = lat, DOY = DOY, f = day_fraction) # W/m2
    # Conversion factors to PAR for direct and diffuse irradiance
    f_dir = waveband_conversion(Itype = :direct,  waveband = :PAR, mode = :power)
    f_dif = waveband_conversion(Itype = :diffuse, waveband = :PAR, mode = :power)
    # Actual irradiance per waveband
    Idir_PAR = f_dir * sky_light[:Idir]
    Idif_PAR = f_dif * sky_light[:Idif]
    # Create the dome of diffuse light
    dome = sky(
        mesh,
        Idir = 0.0, ## No direct solar radiation
        Idif = Idif_PAR, ## Daily Diffuse solar radiation
        nrays_dif = 1_000_000, ## Total number of rays for diffuse solar radiation
        sky_model = StandardSky, ## Angular distribution of solar radiation
        dome_method = equal_solid_angles, # Discretization of the sky dome
        ntheta = 9, ## Number of discretization steps in the zenith angle
        nphi = 12 ## Number of discretization steps in the azimuth angle
    ) 
    # Add direct source
    append!(dome, sky(
            mesh, Idir = Idir_PAR, nrays_dir = 1_000_000,
            Idif = 0.0, nrays_diff = 0, theta_dir = sky_light[:theta], phi_dir = sky_light[:phi]
        )
    )
    
    return dome
end

create_sky (generic function with 1 method)

Create a geometry for the soil surface for light reflection.

function create_soil()
    soil = Rectangle(length = 21.0, width = 21.0)
    rotatey!(soil, π/2) ## To put it in the XY plane
    
    return soil
end

create_soil (generic function with 1 method)

Create a geometry for the entire scene, consisting of our tree and the soil surface. We also add a material for the soil surface.

function create_scene(tree)
    mesh = Mesh(tree)
    soil = create_soil()
    soil_material = Lambertian(τ = 0.0, ρ = 0.21)
    add!(mesh, soil, materials = soil_material)
    
    return mesh
end

create_scene (generic function with 1 method)

Run the ray tracer. This consists of creating the geometry of our scene, creating our light sources, defining the ray tracer with settings of choice and finally running it with trace!.

function run_raytracer!(tree; day_fraction = 0.5)
    mesh = create_scene(tree)
    # the directional light sources created with `SkyDomes.jl` require an "accelerated" mesh to run
    accmesh = accelerate(mesh, acceleration = BVH, rule = SAH{3}(5, 10)) 
    sources = create_sky(day_fraction; mesh = accmesh)
    # note we set `nx` and `ny` to 0: these are grid cloning parameters used to minimize the boundary effects of only simulating a part of a forest. however, as we simulate only a single tree, this is not applicable here.
    settings = RTSettings(pkill = 0.9, maxiter = 4, nx = 0, ny = 0,
                          parallel = true, verbose = false)
    raytracer = RayTracer(mesh, sources; settings)
    trace!(raytracer)
    
    return nothing
end

run_raytracer! (generic function with 1 method)

Running the ray tracer as we solve our differential equations is not computationally feasible. Instead, we run the ray tracer beforehand for a number of times throughout the day and define an interpolation function for every leaf that calculates the incoming PAR for a given time based on these PAR samples.

function precalculate_PAR!(tree; Δf = 0.05)
    # remove potential existing PAR samples
    for node in getnodes(tree)
        if getstructmod(node) == :Leaf
            empty!(data(node).PAR_samples)
        end
    end

    # generate new PAR samples and divide by surface area
    for day_fraction in 0.0:Δf:1.0
        run_raytracer!(tree; day_fraction)
        for node in getnodes(tree)
            if getstructmod(node) == :Leaf
                PAR = data(node).mat |> power |> only
                PAR_flux = PAR / 
                    (surface_area(Cuboid(), data(node).D) * 1e-4) # cm^2 to m^2
                push!(data(node).PAR_samples, PAR_flux)
            end
        end
    end

    # interpolate between samples
    for node in getnodes(tree)
        if getstructmod(node) == :Leaf
            PAR_interpolation = LinearInterpolation(
                data(node).PAR_samples,
                0.0:Δf:1.0,
                extrapolation = ExtrapolationType.Constant # use constant boundary values for extrapolation (here always 0)
            );
            
            PAR_func(t, t_sunrise, t_sunset) = 
                (t % 24 - t_sunrise) / (t_sunset - t_sunrise) |>
                f -> clamp(f, 0.0, 1.0) |>
                PAR_interpolation
            
            data(node).PAR_func = PAR_func
        end
    end
    
    return nothing
end

precalculate_PAR! (generic function with 1 method)

Running the ray tracer

Let's run the ray tracer and visualise the incoming PAR over time for our leaf.

precalculate_PAR!(shoot_graph)

begin
    plot(xlims = (0.0, 24.0), xlabel = "Time (h)", ylabel = "PAR (W/m²)", legend = false)
    for node in getnodes(shoot_graph)
        if getstructmod(node) == :Leaf
            plot!(t -> data(node).PAR_func(t, 8.0, 20.0))
        end
    end
    plot!()
end

Connecting the separate graphs

We finalize the structural definition by connecting all the separate graphs. The graphs are connected as follows:

Shoots and roots: the base of the shoots is connected to the base of the root system.
Roots and soil: each root segment is connected to the soil voxel they are physically inside of, which is done by comparing the coordinates of the root segments and soil voxels. Complicated connections like this can be defined by passing a function to intergraph_connections for the graphs in question that takes a node of each graph and returns whether they should be connected.
Shoots and air: each leaf node is connected to the single air node.
Soil and air: each node of the top soil layer is connected to the air to simulate evaporation.

graphs = [shoot_graph, root_graph, soil_graph, Air()];

function is_connected_root_soil(root, soil)
    soil_center = [soil.x, soil.y, soil.z]
    coord_inside_voxel = [
        soil_center[i]-vs/2 <= data(root).coords[i] < soil_center[i]+vs/2
        for i in eachindex(data(root).coords)
    ]
    are_nodes_connected = all(coord_inside_voxel)
    return are_nodes_connected
end

is_connected_root_soil (generic function with 1 method)

intergraph_connections = [
    (1, 2) => (getnodes(shoot_graph)[1], getnodes(root_graph)[1]),
    (2, 3) => is_connected_root_soil,
    (1, 4) => (:Leaf, :Air),
    (3, 4) => (soil_graph[:, :, 1], :Air)
];

plantstructure = PlantStructure(graphs, intergraph_connections);

plotstructure(plantstructure)

Functional definition

For the definition of the functional processes, we return to PlantModules.jl territory. Our structural modules are assigned the same classic sets of functional modules as in previous tutorials, with two exceptions:

We don't simulate hydraulics for the meristem, which is only a thin layer of cells.
We simulate light interception and carbon assimilation in the leaves, with a guest appearance of the PlantBioPhysics.jl package (also from the VirtualPlantLab ecosystem) for the calculation of carbon assimilation.

Functional modules

import PlantModules: t, d

Our meristem is part of the graph, so it needs to respect the PlantModules.multi_connection_eqs function, which defines how nodes interact with their neighbours. The default value for this function simply states that the net water inflow \(ΣF\) of a node equals the sum of the water flows \(F\) between the node and each of its neighbors. Therefore, if we want the meristem to be hydraulically inactive, we need to set both its water inflow \(ΣF\) and each of its water flows \(F\) to \(0\).

function inactive_hydraulic_module(; name)
    @variables ΣF(t)
    eqs = [d(ΣF) ~ 0]

    return System(eqs, t; name)
end

inactive_hydraulic_module (generic function with 1 method)

function inactive_hydraulic_connection(; name)
    @variables F(t) = 0.0

    eqs = [d(F) ~ 0]
    # `get_connection_eqset` defines how variables of the connection relate to variables of two nodes it connects, but as we use no variables here that correspond with those from the nodes, the set is simply empty 
    get_connection_eqset(node_MTK, nb_node_MTK, connection_MTK) = []
    
    return System(eqs, t; name), get_connection_eqset
end

inactive_hydraulic_connection (generic function with 1 method)

Based on our ray tracer, we can calculate the incoming PAR for every leaf at any time of day. However, we still need to model how incoming PAR relates to carbon assimilation. We will use the established carbon assimilation and stomatal conducatance models available in PlantBiophysics.jl for this.

function get_assimilation_rate(PAR_flux, T)
    meteo = Atmosphere(
        T = T, Wind = 1.0, P = 101.3, Rh = 0.65
    )
    m = ModelList(
        Fvcb(), # calculate CO2 assimilation rate
        Medlyn(0.03, 0.92), # calculate stomatal conductance, see https://onlinelibrary.wiley.com/doi/epdf/10.1111/j.1365-2486.2010.02375.x
        status = (Tₗ = meteo[:T], Cₛ = meteo[:Cₐ], Dₗ = meteo[:VPD], RI_PAR_f = meteo[:Ri_PAR_f], aPPFD = PAR_flux)
    )
    run!(m, meteo)
    assimilation_rate = m[:A][1] # extract result of the first (and only) timestep
    assimilation_rate = max(0, assimilation_rate) # set negative values to 0
    return assimilation_rate
end

get_assimilation_rate (generic function with 1 method)

When we want to use more complex functions such as this inside of our functional modules, it's often a good idea to reduce their computational cost if possible. For this example, we can simply replace the function with another linear interpolation, where we set the temperature to the value we will use in our simulations.

interpolation_range = 0:1000

0:1000

get_assimilation_rate_interpolation = LinearInterpolation(
    get_assimilation_rate.(interpolation_range, 25.0),
    interpolation_range,
    extrapolation = ExtrapolationType.Extension # smoothly extend interpolation for extrapolation
);

Now we can define the photosynthesis module itself. It is rather straightforward, as we have already defined our functions for getting incoming PAR based on the time of day and carbon assimilation based on the incoming PAR. Note that to access the PAR_func defined in the graph nodes we treat PAR_func as any other parameter: we include it in the inputs of our functional module and assign it a default parameter value later on. This way, our functional module will grab the value for PAR_func defined in the graph nodes during model creation.

function photosynthesis_module(; name, M, shape, M_c, PAR_func, t_sunrise, t_sunset)
    @constants (
        uc = (10^-6 * 10^-4 * 60^2), [description = "Unit conversion from (µmol / m^2 / s) to (mol / cm^2 / hr)"],
            # the output from PlantBiophysics.jl is in different units than we use for our ODEs, so we need to change this
    )
    @parameters (
        M_c = M_c, [description = "Rate of carbon consumption"], # hr^-1
    )
    @variables (
        M(t) = M, [description = "Osmotically active metabolite content"], 
            # mol / cm^3
        PF(t), [description = "Incoming PAR flux"], # W / m^2
        A(t), [description = "Carbon assimilation rate"], # mol / cm^2 / hr
        A_V(t), [description = "Volumetric carbon assimilation rate"], # mol / cm^3 / hr
        D(t)[1:getdimensionality(shape)], [description = "Dimensions of compartment"], # cm
    )

    eqs = [
        PF ~ PAR_func(t, t_sunrise, t_sunset)
        A ~ uc * get_assimilation_rate_interpolation(PF)
        A_V ~ A * surface_area(shape, D) / PlantModules.volume(shape, D)
        d(M) ~ A_V - M_c*M
    ]
    return System(eqs, t; name, checks = false)
end

photosynthesis_module (generic function with 1 method)

Coupling

We couple structural and functional modules as discussed at the start of this section.

module_coupling = Dict(
    :Meristem => [inactive_hydraulic_module],
    :Bud => [hydraulic_module, constant_carbon_module, K_module],
    :Node => [hydraulic_module, constant_carbon_module, K_module],
    :BudNode => [hydraulic_module, constant_carbon_module, K_module],
    :Internode => [hydraulic_module, constant_carbon_module, K_module],
    :Leaf => [hydraulic_module, photosynthesis_module, K_module],

    :Root => [hydraulic_module, constant_carbon_module, K_module],
    
    :Soil => [environmental_module, Ψ_soil_module, constant_K_module],
    :Air => [environmental_module, Ψ_air_module, constant_K_module],
);

connecting_modules = Dict(
    (:Soil, :Soil) => constant_hydraulic_connection,
    (:Soil, :Root) => constant_hydraulic_connection,
    (:Root, :Root) => hydraulic_connection,
    (:Root, :Internode) => hydraulic_connection,
    (:Internode, :Node) => hydraulic_connection,
    (:Node, :Bud) => hydraulic_connection,
    (:Node, :BudNode) => hydraulic_connection,
    (:BudNode, :Internode) => hydraulic_connection,
    (:Node, :Leaf) => hydraulic_connection,
    (:Internode, :Meristem) => inactive_hydraulic_connection,
    (:Leaf, :Air) => daynight_hydraulic_connection,
    (:Soil, :Air) => constant_hydraulic_connection
);

plantcoupling = PlantCoupling(; module_coupling, connecting_modules);

Parameters

Parameter specification follows the usual steps of assigning the correct shapes and initial dimensions to the plant parts and setting the water capacity W_max of our soil, initial relative water content W_r of the air, and hydraulic conductivity K of the air. We also lower the hydraulic conductivity between soil and air to a more realistic value for direct evaporation from the soil.

default_changes = Dict(
    :PAR_func => (x...) -> error("A leaf node exists with no defined PAR function.")
)

Dict{Symbol, Main.var"#56#57"{typeof(error)}} with 1 entry:
  :PAR_func => #56

module_defaults = Dict(
    :Bud => Dict(:shape => PlantModules.Sphere()),
    :Leaf => Dict(:shape => PlantModules.Cuboid()),
    :Node => Dict(:D => [0.5, 0.1]),
    :BudNode => Dict(:D => [0.5, 0.1]),
    :Soil => Dict(:W_max => 1e4),
    :Air => Dict(:W_r => 0.6, :K => 1e-3)
);

connection_values = Dict(
    (:Soil, :Air) => Dict(:K => 1e-2)
)

Dict{Tuple{Symbol, Symbol}, Dict{Symbol, Float64}} with 1 entry:
  (:Soil, :Air) => Dict(:K=>0.01)

plantparams = PlantParameters(; default_changes, module_defaults, connection_values);

Running the model

Finally, we generate and run the system. We will start with an initial run to generate our solution for the first day as in our previous two tutorials. For the following days, we define a function that performs the following steps:

Update the structures of the shoot and root graphs, the former of which depends on the functional status of the plant evaluated at the end of the previous day.
Run our ray tracer again to recalculate the incoming PAR for our updated shoot structure.
Create an updated PlantStructure using the updated shoot and root graphs.
Add the solution of the previous day's simulation to PlantParameters. This makes it use the final values of all variables as the initial values for the simulation of the next day.
Generate and solve the system.
Store our solution and plantstructure for visualisation of the results.
Based on the day's simulation, update the water contents of bud nodes and the dimensions of the internodes and leaves in our shoot graph. These are no longer used directly by our package, as the initial values they specify are now overwritten by the final variable values specified in the previous day's solution. However, they are used in the next day's rewriting steps and ray tracing, so it is important to update them.

Initial run

system = generate_system(plantstructure, plantcoupling, plantparams);

tspan = (0.0, 24.0);

prob = ODEProblem(system, [], tspan, sparse = true);

sol = solve(prob, FBDF());

Subsequent runs

function run_timestep!(sols, plantstructures, shoot_graph, root_graph)
    # rewrite structures
    rewrite!(shoot_graph)
    rewrite!(root_graph)

    # run ray tracer for new leaves
    precalculate_PAR!(shoot_graph)

    # connect new structures
    graphs = [shoot_graph, root_graph, soil_graph, Air()]
    intergraph_connections = [
        (1, 2) => (getnodes(shoot_graph)[1], getnodes(root_graph)[1]),
        (2, 3) => is_connected_root_soil,
        (1, 4) => (:Leaf, :Air),
        (3, 4) => (soil_graph[:, :, 1], :Air)
    ]
    plantstructure_new = PlantStructure(graphs, intergraph_connections);

    # add solution to parameters - changes initial values to previous solution's end values
    plantparams_new = PlantParameters(; default_changes, module_defaults,
                                      connection_values, sol = sols[end]);

    # generate system and solve
    system_new = generate_system(plantstructure_new, plantcoupling, plantparams_new);
    prob_new = ODEProblem(system_new, [], tspan, sparse = true);
    sol_new = solve(prob_new, FBDF());

    # store data for plotting
    push!(sols, sol_new)
    push!(plantstructures, plantstructure_new)

    # update node data used in rewriting rules, plotting and ray tracing
    for node in getnodes(shoot_graph)
        if getstructmod(node) == :Bud
            var = get_subsystem_variables(
                system_new, shoot_graph, 1, :W, node
            )
            getattributes(node)[:W] = sol_new[var][end]
        elseif getstructmod(node) in [:Leaf, :Internode]
            var = get_subsystem_variables(
                system_new, shoot_graph, 1, :D, node
            )
            getattributes(node)[:D] = sol_new[var][end]
        end
    end
    
    return nothing
end

run_timestep! (generic function with 1 method)

days = 5

sols = [sol];

plantstructures = [plantstructure];

for day in 2:days
    @info "Simulating day $day"
    run_timestep!(sols, plantstructures, shoot_graph, root_graph)
end

Let's quickly inspect the structure of our plant at the end of the simulation period.

render(Mesh(shoot_graph))

plotstructure(root_graph)

plotstructure(plantstructures[end])

Results

Finally, we can visualize the results of our simulation over the entire time period. We will make the following plots to show the effects of the new functionalities discussed in this tutorial

Bud water content over time, illustrating the structural growth of our plant as new buds grow and old buds grow into branches.
Soil water content over time, illustrating the effect of discretizing the soil into multiple compartments.
The leaf assimilation rate over time, illustrating the different trends of incoming PAR for different leaves.

Plotting has become technically more complex now that we have multiple plantstructures and corresponding solutions, but we can simply pass them to plotgraph as vectors and use the plotting function as per usual. Plotting soil water content over time is more complicated, however, because we want to categorize them based on whether they are directly below the plant.