Expanding on the multiscale simulation

Let's build on the previous example and add some other organ growth, as well as some very mild coupling between the two.

You can find the full script for this simulation in the ToyMultiScalePlantModel subfolder of the examples folder.

• Expanding on the multiscale simulation
• • Setup
  • Adding roots to our plant
  • Root models
  • • Water absorption
    • Root growth
  • Updating other models to account for water
  • • Resource storage
    • Internode creation
  • Updating the mapping
  • Running the simulation

Setup

Once again, with a properly set-up Julia environment:

using PlantSimEngine
using PlantSimEngine.Examples
using PlantMeteo
using MultiScaleTreeGraph
using CSV, DataFrames

PlantSimEngine.@process "leaf_carbon_capture" verbose = false

struct ToyLeafCarbonCaptureModel<: AbstractLeaf_Carbon_CaptureModel end

function PlantSimEngine.inputs_(::ToyLeafCarbonCaptureModel)
    NamedTuple() # No inputs
end

function PlantSimEngine.outputs_(::ToyLeafCarbonCaptureModel)
    (carbon_captured=0.0,)
end

function PlantSimEngine.run!(::ToyLeafCarbonCaptureModel, models, status, meteo, constants, extra)
    status.carbon_captured = 40
end

function get_n_leaves(node::MultiScaleTreeGraph.Node)
    root = MultiScaleTreeGraph.get_root(node)
    nleaves = length(MultiScaleTreeGraph.traverse(root, x->1, symbol="Leaf"))
    return nleaves
end

get_n_leaves (generic function with 1 method)

Adding roots to our plant

We'll add a root that extracts water and adds it to the stock. Initial water stocks are low, so root growth is prioritized, then the plant also grows leaves and a new internode like it did before. Roots only grow up to a certain point, and don't branch.

This leads to adding a new scale, "Root" to the mapping, as well as two more models, one for water absorption, the other for root growth. Other models are updated here and there to account for water. The carbon capture model remains unchanged, and so is the get_n_leaves helper function.

Root models

Water absorption

Let's implement a very fake model of root water absorption. It'll capture the amount of precipitation in the weather data multiplied by some assimilation factor.

PlantSimEngine.@process "water_absorption" verbose = false

struct ToyWaterAbsorptionModel <: AbstractWater_AbsorptionModel
end

PlantSimEngine.inputs_(::ToyWaterAbsorptionModel) = (root_water_assimilation=1.0,)
PlantSimEngine.outputs_(::ToyWaterAbsorptionModel) = (water_absorbed=0.0,)

function PlantSimEngine.run!(m::ToyWaterAbsorptionModel, models, status, meteo, constants=nothing, extra=nothing)
    status.water_absorbed = meteo.Precipitations * status.root_water_assimilation
end

Root growth

The root growth model is similar to the internode growth one : it checks for a water threshold and that there is enough carbon, and adds a new organ to the MTG if the maximum length hasn't been reached.

It also makes use of a couple of helper functions to find the end root and compute root length :

function get_root_end_node(node::MultiScaleTreeGraph.Node)
    root = MultiScaleTreeGraph.get_root(node)
    return MultiScaleTreeGraph.traverse(root, x->x, symbol="Root", filter_fun = MultiScaleTreeGraph.isleaf)
end

function get_roots_count(node::MultiScaleTreeGraph.Node)
    root = MultiScaleTreeGraph.get_root(node)
    return length(MultiScaleTreeGraph.traverse(root, x->x, symbol="Root"))
end

PlantSimEngine.@process "root_growth" verbose = false

struct ToyRootGrowthModel{T} <: AbstractRoot_GrowthModel
    water_threshold::T
    carbon_root_creation_cost::T
    root_max_len::Int
end

PlantSimEngine.inputs_(::ToyRootGrowthModel) = (water_stock=0.0,carbon_stock=0.0,)
PlantSimEngine.outputs_(::ToyRootGrowthModel) = (carbon_root_creation_consumed=0.0,)

function PlantSimEngine.run!(m::ToyRootGrowthModel, models, status, meteo, constants=nothing, extra=nothing)
    if status.water_stock < m.water_threshold && status.carbon_stock > m.carbon_root_creation_cost

        root_end = get_root_end_node(status.node)

        if length(root_end) != 1
            throw(AssertionError("Couldn't find MTG leaf node with symbol \"Root\""))
        end
        root_len = get_roots_count(root_end[1])
        if root_len < m.root_max_len
            st = add_organ!(root_end[1], extra, "<", "Root", 2, index=1)
            status.carbon_root_creation_consumed = m.carbon_root_creation_cost
        end
    else
        status.carbon_root_creation_consumed = 0.0
    end
end

Updating other models to account for water

Resource storage

Water absorbed must now be accumulated, and root carbon creation costs taken into account.

PlantSimEngine.@process "resource_stock_computation" verbose = false

struct ToyStockComputationModel <: AbstractResource_Stock_ComputationModel
end

PlantSimEngine.inputs_(::ToyStockComputationModel) =
(water_absorbed=0.0,carbon_captured=0.0,carbon_organ_creation_consumed=0.0,carbon_root_creation_consumed=0.0)

PlantSimEngine.outputs_(::ToyStockComputationModel) = (water_stock=-Inf,carbon_stock=-Inf)

function PlantSimEngine.run!(m::ToyStockComputationModel, models, status, meteo, constants=nothing, extra=nothing)
    status.water_stock += sum(status.water_absorbed)
    status.carbon_stock += sum(status.carbon_captured) - sum(status.carbon_organ_creation_consumed) - sum(status.carbon_root_creation_consumed)
end

Internode creation

The minor change is that new organs are now created only if the water stock is above a given threshold.

struct ToyCustomInternodeEmergence{T} <: AbstractOrgan_EmergenceModel
    TT_emergence::T
    carbon_internode_creation_cost::T
    leaf_surface_area::T
    leaves_max_surface_area::T
    water_leaf_threshold::T
end

ToyCustomInternodeEmergence(;TT_emergence=300.0, carbon_internode_creation_cost=200.0, leaf_surface_area=3.0,leaves_max_surface_area=100.0,
water_leaf_threshold=30.0) = ToyCustomInternodeEmergence(TT_emergence, carbon_internode_creation_cost, leaf_surface_area, leaves_max_surface_area, water_leaf_threshold)

PlantSimEngine.inputs_(m::ToyCustomInternodeEmergence) = (TT_cu=0.0,water_stock=0.0, carbon_stock=0.0)
PlantSimEngine.outputs_(m::ToyCustomInternodeEmergence) = (TT_cu_emergence=0.0, carbon_organ_creation_consumed=0.0)

function PlantSimEngine.run!(m::ToyCustomInternodeEmergence, models, status, meteo, constants=nothing, sim_object=nothing)

    leaves_surface_area = m.leaf_surface_area * get_n_leaves(status.node)
    status.carbon_organ_creation_consumed = 0.0

    if leaves_surface_area > m.leaves_max_surface_area
        return nothing
    end

    # if water levels are low, prioritise roots
    if status.water_stock < m.water_leaf_threshold
        return nothing
    end

    # if not enough carbon, no organ creation
    if status.carbon_stock < m.carbon_internode_creation_cost
        return nothing
    end

    if length(MultiScaleTreeGraph.children(status.node)) == 2 &&
        status.TT_cu - status.TT_cu_emergence >= m.TT_emergence
        status_new_internode = add_organ!(status.node, sim_object, "<", "Internode", 2, index=1)
        add_organ!(status_new_internode.node, sim_object, "+", "Leaf", 2, index=1)
        add_organ!(status_new_internode.node, sim_object, "+", "Leaf", 2, index=1)

        status_new_internode.TT_cu_emergence = m.TT_emergence - status.TT_cu
        status.carbon_organ_creation_consumed = m.carbon_internode_creation_cost
    end

    return nothing
end

Updating the mapping

The resource storage and internode emergence models now need a couple of extra water-related mapped variables. The "Root" organ is added to the mapping with its own models. New parameters need to be initialized.

mapping = Dict(
"Scene" => ToyDegreeDaysCumulModel(),
"Plant" => (
    MultiScaleModel(
        model=ToyStockComputationModel(),
        mapped_variables=[
            :carbon_captured=>["Leaf"],
            :water_absorbed=>["Root"],
            :carbon_root_creation_consumed=>["Root"],
            :carbon_organ_creation_consumed=>["Internode"]

        ],
        ),
        Status(water_stock = 0.0, carbon_stock = 0.0)
    ),
"Internode" => (
        MultiScaleModel(
            model=ToyCustomInternodeEmergence(),#TT_emergence=20.0),
            mapped_variables=[:TT_cu => "Scene",
            PreviousTimeStep(:water_stock)=>"Plant",
            PreviousTimeStep(:carbon_stock)=>"Plant"],
        ),
        Status(carbon_organ_creation_consumed=0.0),
    ),
"Root" => ( MultiScaleModel(
            model=ToyRootGrowthModel(10.0, 50.0, 10),
            mapped_variables=[PreviousTimeStep(:carbon_stock)=>"Plant",
            PreviousTimeStep(:water_stock)=>"Plant"],
        ),
            ToyWaterAbsorptionModel(),
            Status(carbon_root_creation_consumed=0.0, root_water_assimilation=1.0),
            ),
"Leaf" => ( ToyLeafCarbonCaptureModel(),),
)

Dict{String, Any} with 5 entries:
  "Internode" => (MultiScaleModel{ToyCustomInternodeEmergence{Float64}, Vector{…
  "Root"      => (MultiScaleModel{ToyRootGrowthModel{Float64}, Vector{Pair{Unio…
  "Scene"     => ToyDegreeDaysCumulModel(0.0, 10.0, 43.0)
  "Plant"     => (MultiScaleModel{ToyStockComputationModel, Vector{Pair{Union{S…
  "Leaf"      => (ToyLeafCarbonCaptureModel(),)

Running the simulation

Running this new simulation is almost the same as before. The weather data is unchanged, but a new "Root" node was added to the MTG.

mtg = MultiScaleTreeGraph.Node(MultiScaleTreeGraph.NodeMTG("/", "Scene", 1, 0))
    plant = MultiScaleTreeGraph.Node(mtg, MultiScaleTreeGraph.NodeMTG("+", "Plant", 1, 1))

    internode1 = MultiScaleTreeGraph.Node(plant, MultiScaleTreeGraph.NodeMTG("/", "Internode", 1, 2))
    MultiScaleTreeGraph.Node(internode1, MultiScaleTreeGraph.NodeMTG("+", "Leaf", 1, 2))
    MultiScaleTreeGraph.Node(internode1, MultiScaleTreeGraph.NodeMTG("+", "Leaf", 1, 2))

    internode2 = MultiScaleTreeGraph.Node(internode1, MultiScaleTreeGraph.NodeMTG("<", "Internode", 1, 2))
    MultiScaleTreeGraph.Node(internode2, MultiScaleTreeGraph.NodeMTG("+", "Leaf", 1, 2))
    MultiScaleTreeGraph.Node(internode2, MultiScaleTreeGraph.NodeMTG("+", "Leaf", 1, 2))

    plant_root_start = MultiScaleTreeGraph.Node(
        plant,
        MultiScaleTreeGraph.NodeMTG("+", "Root", 1, 3),
    )

meteo_day = CSV.read(joinpath(pkgdir(PlantSimEngine), "examples/meteo_day.csv"), DataFrame, header=18)

outs = run!(mtg, mapping, meteo_day)
mtg

/ 1: Scene
└─ + 2: Plant
   ├─ / 3: Internode
   │  ├─ + 4: Leaf
   │  ├─ + 5: Leaf
   │  └─ < 6: Internode
   │     ├─ + 7: Leaf
   │     ├─ + 8: Leaf
   │     └─ < 19: Internode
   │        ├─ + 20: Leaf
   │        ├─ + 21: Leaf
   │        └─ < 22: Internode
   │           ├─ + 23: Leaf
   │           ├─ + 24: Leaf
   │           └─ < 25: Internode
   │              ├─ + 26: Leaf
   │              ├─ + 27: Leaf
   │              └─ < 28: Internode
   │                 ├─ + 29: Leaf
   │                 ├─ + 30: Leaf
   │                 └─ < 31: Internode
   │                    ├─ + 32: Leaf
   │                    ├─ + 33: Leaf
   │                    └─ < 34: Internode
…

And that's it !

...Or is it ?

If you inspect the code and output data closely, you may notice some distinctive problems with the way the simulation runs... Some things aren't quite right. If you wish to know more, onwards to the next chapter: Fixing bugs in the plant simulation