GitHub

Converting a single-scale simulation to multi-scale

A single-scale simulation can be turned into a 'pseudo-multi-scale' simulation by providing a simple multi-scale tree graph, and declaring a mapping linking all models to a unique scale level.

This page showcases how to do the conversion, and then adds a model at a new scale to make the simulation genuinely multi-scale.

The full script for the example can be found in the examples folder, here

Converting the ModelList to a multi-scale mapping

For example, let's return to the ModelList coupling a light interception model, a Leaf Area Index model, and a carbon biomass increment model that was discussed in the Model switching subsection: 

using PlantMeteo
using PlantSimEngine
using PlantSimEngine.Examples
using CSV

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

models_singlescale = ModelList(
    ToyLAIModel(),
    Beer(0.5),
    ToyRUEGrowthModel(0.2),
    status=(TT_cu=cumsum(meteo_day.TT),),
)

outputs_singlescale = run!(models_singlescale, meteo_day)
TimeStepTable{Status{(:TT_cu, :LAI, :aPPFD,...}(365 x 5):
╭─────┬──────────┬────────────┬───────────┬────────────┬───────────────────╮
│ Row │    TT_cu │        LAI │     aPPFD │    biomass │ biomass_increment │
│     │  Float64 │    Float64 │   Float64 │    Float64 │           Float64 │
├─────┼──────────┼────────────┼───────────┼────────────┼───────────────────┤
│   1 │      0.0 │ 0.00554988 │ 0.0396961 │ 0.00793922 │        0.00793922 │
│   2 │      0.0 │ 0.00554988 │   0.02173 │  0.0122852 │          0.004346 │
│   3 │      0.0 │ 0.00554988 │ 0.0314899 │  0.0185832 │        0.00629798 │
│   4 │      0.0 │ 0.00554988 │ 0.0390834 │  0.0263999 │        0.00781668 │
│   5 │      0.0 │ 0.00554988 │ 0.0454514 │  0.0354902 │        0.00909028 │
│   6 │      0.0 │ 0.00554988 │ 0.0472677 │  0.0449437 │        0.00945354 │
│   7 │      0.0 │ 0.00554988 │   0.04346 │  0.0536357 │        0.00869201 │
│   8 │      0.0 │ 0.00554988 │ 0.0469832 │  0.0630324 │        0.00939665 │
│   9 │      0.0 │ 0.00554988 │ 0.0291703 │  0.0688664 │        0.00583406 │
│  10 │      0.0 │ 0.00554988 │ 0.0140052 │  0.0716675 │        0.00280105 │
│  11 │      0.0 │ 0.00554988 │ 0.0505283 │  0.0817731 │         0.0101057 │
│  12 │      0.0 │ 0.00554988 │ 0.0405277 │  0.0898787 │        0.00810554 │
│  13 │   0.5625 │ 0.00557831 │ 0.0297814 │   0.095835 │        0.00595629 │
│  14 │ 0.945833 │ 0.00559777 │ 0.0433269 │     0.1045 │        0.00866538 │
│  15 │ 0.979167 │ 0.00559946 │ 0.0470271 │   0.113906 │        0.00940542 │
│  ⋮  │    ⋮     │     ⋮      │     ⋮     │     ⋮      │         ⋮         │
╰─────┴──────────┴────────────┴───────────┴────────────┴───────────────────╯
                                                            350 rows omitted

Those models all operate on a simplified model of a single plant, without any organ-local information. We can therefore consider them to be working at the 'whole plant' scale. Their variables also operate at that "plant" scale, so there is no need to map any variable to other scales.

We can therefore convert this into the following mapping : 

mapping = Dict(
"Plant" => (
   ToyLAIModel(),
    Beer(0.5),
    ToyRUEGrowthModel(0.2),
    Status(TT_cu=cumsum(meteo_day.TT),)
    ),
)
Dict{String, Tuple{PlantSimEngine.Examples.ToyLAIModel, PlantSimEngine.Examples.Beer{Float64}, PlantSimEngine.Examples.ToyRUEGrowthModel{Float64}, Status{(:TT_cu,), Tuple{Base.RefValue{Vector{Float64}}}}}} with 1 entry:
  "Plant" => (ToyLAIModel(8.0, 800, 110.0, 1500, 20.0), Beer{Float64}(0.5), Toy…

Note the slight difference in syntax for the Status. This is due to an implementation quirk (sorry).

Adding a new package for our plant graph

None of these models operate on a multi-scale tree graph, either. There is no concept of organ creation or growth. We still need to provide a multi-scale tree graph to a multi-scale simulation, so we can -for now- declare a very simple MTG, with a single node:

using MultiScaleTreeGraph

mtg = MultiScaleTreeGraph.Node(MultiScaleTreeGraph.NodeMTG("/", "Plant", 0, 0),)
/ 1: Plant
Note

You will need to add the MultiScaleTreeGraph package to your environment. See Installing and running PlantSimEngine if you are not yet comfortable with Julia or need a refresher.

Running the multi-scale simulation ?

We now have almost everything we need to run the multiscale simulation.

This first conversion step can be a starting point for a more elaborate multi-scale simulation.

The signature of the run! function in multi-scale differs slightly from the ModelList version : 

out_multiscale = run!(mtg, mapping, meteo_day)

(Some of the optional arguments also change slightly)

Unfortunately, there is one caveat. Passing in a vector through the Status field is still possible in multi-scale mode, but requires a little more advanced tinkering with the mapping, as it generates a custom model under the hood and the implementation is experimental and less user-friendly.

If you are keen on going down that path, you can find a detailed example here, but we don't recommend it for beginners.

What we'll do instead, is write our own model provide the thermal time per timestep as a variable, instead of as a single vector in the Status.

Our 'pseudo-multiscale' first approach will therefore turn into a genuine multi-scale simulation.

Adding a second scale

Let's have a model provide the Cumulated Thermal Time to our Leaf Area Index model, instead of initializing it through the Status.

Let's instead implement our own ToyTT_cuModel.

TT_cu model implementation

This model doesn't require any outside data or input variables, it only operates on the weather data and outputs our desired TT_cu. The implementation doesn't require any advanced coupling and is very straightforward.

PlantSimEngine.@process "tt_cu" verbose = false

struct ToyTt_CuModel <: AbstractTt_CuModel
end

function PlantSimEngine.run!(::ToyTt_CuModel, models, status, meteo, constants, extra=nothing)
    status.TT_cu +=
        meteo.TT
end

function PlantSimEngine.inputs_(::ToyTt_CuModel)
    NamedTuple() # No input variables
end

function PlantSimEngine.outputs_(::ToyTt_CuModel)
    (TT_cu=0.0,)
end
Note

The only accessible variables in the run! function via the status are the ones that are local to the "Scene" scale. This isn't explicit at first glance, but very important to keep in mind when developing models, or using them at different scales. If variables from other scales are required, then they need to be mapped via a MultiScaleModel, or sometimes a more complex coupling is necessary.

Linking the new TT_cu model to a scale in the mapping

We now have our model implementation. How does it fit into our mapping ?

Our new model doesn't really relate to a specific organ of our plant. In fact, this model doesn't represent a physiological process of the plant, but rather an environmental process affecting its physiology. We could therefore have it operate at a different scale unrelated to the plant, which we'll call "Scene". This makes sense.

Note that we now need to add a "Scene" node to our Multi-scale Tree Graph, otherwise our model will not run, since no other model calls it and "Plant" nodes will only call models at the "Plant" scale. See Empty status vectors in multi-scale simulations for more details.

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

Mapping between scales : the MultiScaleModel wrapper

The cumulated thermal time (:TT_cu) which was previously provided to the LAI model as a simulation parameter now needs to be mapped from the "Scene" scale level.

This is done by wrapping our ToyLAIModel in a dedicated structure called a MultiScaleModel. A MultiScaleModel requires two keyword arguments : model, indicating the model for which some variables are mapped, and mapped_variables, indicating which scale link to which variables, and potentially renaming them.

There can be different kinds of variable mapping with slightly different syntax, but in our case, only a single scalar value of the TT_cu is passed from the "Scene" to the "Plant" scale.

This gives us the following declaration with the MultiScaleModel wrapper for our LAI model: 

MultiScaleModel(
            model=ToyLAIModel(),
            mapped_variables=[
                :TT_cu => "Scene",
            ],
        )
MultiScaleModel{PlantSimEngine.Examples.ToyLAIModel, Vector{Pair{Union{Symbol, PreviousTimeStep}, Union{Pair{String, Symbol}, Vector{Pair{String, Symbol}}}}}}(PlantSimEngine.Examples.ToyLAIModel(8.0, 800, 110.0, 1500, 20.0), Pair{Union{Symbol, PreviousTimeStep}, Union{Pair{String, Symbol}, Vector{Pair{String, Symbol}}}}[:TT_cu => ("Scene" => :TT_cu)])

and the new mapping with two scales:

mapping_multiscale = Dict(
    "Scene" => ToyTt_CuModel(),
    "Plant" => (
        MultiScaleModel(
            model=ToyLAIModel(),
            mapped_variables=[
                :TT_cu => "Scene",
            ],
        ),
        Beer(0.5),
        ToyRUEGrowthModel(0.2),
    ),
)
Dict{String, Any} with 2 entries:
  "Scene" => ToyTt_CuModel()
  "Plant" => (MultiScaleModel{ToyLAIModel, Vector{Pair{Union{Symbol, PreviousTi…

Running the multi-scale simulation

We can then run the multiscale simulation, with our two-node MTG :

outputs_multiscale = run!(mtg_multiscale, mapping_multiscale, meteo_day)
Dict{String, Dict{Symbol, Vector}} with 2 entries:
  "Scene" => Dict(:TT_cu=>[[0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.…
  "Plant" => Dict(:TT_cu=>[[0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.…

Comparing outputs between single- and multi-scale

The outputs structures are slightly different : multi-scale outputs are indexed by scale, and a variable has a value for every node of the scale it operates at (for instance, there would be a "leaf_surface" value for every leaf in a plant), stored in an array.

In our simple example, we only have one MTG scene node and one plant node, so the arrays for each variable in the multi-scale output only contain one value.

We can access the output variables at the "Scene" scale by indexing our outputs:

outputs_multiscale["Scene"]
Dict{Symbol, Vector} with 2 entries:
  :TT_cu => [[0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.…
  :node  => Vector{Node{NodeMTG, Dict{Symbol, Any}}}[[/ 1: Scene…

and then the computed :TT_cu:

outputs_multiscale["Scene"][:TT_cu]
365-element Vector{Vector{Float64}}:
 [0.0]
 [0.0]
 [0.0]
 [0.0]
 [0.0]
 [0.0]
 [0.0]
 [0.0]
 [0.0]
 [0.0]
 ⋮
 [2177.904166666665]
 [2178.4249999999984]
 [2179.729166666665]
 [2180.4749999999985]
 [2181.1458333333317]
 [2183.941666666665]
 [2189.2666666666646]
 [2192.895833333331]
 [2193.8166666666643]

As you can see, it is a Vector{Vector{T}}, whereas our single-scale output is a Vector{T}:

outputs_singlescale.TT_cu
365-element Vector{Float64}:
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    ⋮
 2177.9041666666662
 2178.4249999999997
 2179.729166666666
 2180.4749999999995
 2181.145833333333
 2183.941666666666
 2189.2666666666664
 2192.895833333333
 2193.816666666666

To compare them value-by-value, we can flatten the multiscale Vector and then do a piecewise approximate equality test :

computed_TT_cu_multiscale = collect(Base.Iterators.flatten(outputs_multiscale["Scene"][:TT_cu]))

for i in 1:length(computed_TT_cu_multiscale)
    if !(computed_TT_cu_multiscale[i] ≈ outputs_singlescale.TT_cu[i])
        println(i)
    end
end

or equivalently, with broadcasting, we can write :

is_approx_equal = length(unique(computed_TT_cu_multiscale .≈ outputs_singlescale.TT_cu)) == 1
true
Note

You may be wondering why we check for approximate equality rather than strict equality. The reason for that is due to floating-point accumulation errors, which are discussed in more detail in Floating-point considerations.

ToyDegreeDaysCumulModel

There is a model able to provide Thermal Time based on weather temperature data, ToyDegreeDaysCumulModel, which can also be found in the examples folder.

We didn't make use of it here for learning purposes. It also computes a thermal time based on default parameters that don't correspond to the thermal time in the example weather data, so results differ from the thermal time already present in the weather data without tinkering with the parameters.