LULUC implementation

[marcadella]

• Massive update to Biomee implementation using OOP and linked list to store cohorts.
• Implementation LULUC for Biomee #272
• LULUC vignette #273
• Correct state transitions to match states #274
• Biome LUH2 interface #286

[fabern]

Using colMeans for the test averages over time, right?
This is much less strict than all daily/yearly values to be correct.
I suggest to remove the colMeans to have stricter testing.
Has this been tried out previously and was rejected? If so for what reasons?

[fabern]

Furthermore, instead of testthat::expect_true testthat::expect_equal() yields more precise error messages, indicating which column (e.g. GPP, or Tree density) disagrees. I argue that this speeds up debugging since this info will already be available from the output of the CI run.

See below the error messages resulting from these two functions:

df_reference <- data.frame(LE=c(1,2),     GPP=c(2,3), Density_treeha=c(2,3))
df_to_test   <- data.frame(LE=c(1.001,2), GPP=c(20.001,30.001), Density_treeha=c(2,3))

testthat::expect_equal(
  tolerance = 1e-4,
  df_to_test,
  df_reference)
#> Error: `df_to_test` not equal to `df_reference`.
#> Component "LE": Mean relative difference: 0.000999001
#> Component "GPP": Mean relative difference: 0.900004

testthat::expect_equal(
  tolerance = 0.01,
  df_to_test,
  df_reference)
#> Error: `df_to_test` not equal to `df_reference`.
#> Component "GPP": Mean relative difference: 0.900004

testthat::expect_true(
  all.equal(
    colMeans(df_to_test), 
    colMeans(df_reference),
    tolerance = 1e-4))
#> Error: all.equal(colMeans(df_to_test), colMeans(df_reference), tolerance = 1e-04) is not TRUE
#> 
#> `actual` is a character vector ('Mean relative difference: 0.8490651')
#> `expected` is a logical vector (TRUE)

^{Created on 2025-02-27 with reprex v2.1.0}

[fabern]

Illustrating it with the real output, and mocked up errors:

library(rsofun)
library(dplyr)
library(tidyr)
df_reference2 <- rsofun::biomee_p_model_output$data[[1]]$output_daily_tile |> tidyr::tibble()
df_to_test2   <- df_reference2 |>
  dplyr::mutate(Evap     = Evap * (1 + runif(n(), min = -1e-3, max = 1e-3)),
                GPP      = GPP  * (1 + runif(n(), min = -1e-4, max = 1e-4)),
                slowSOM  = slowSOM  * (1 + runif(n(), min = -1e-4, max = 1e-4)))

testthat::expect_true(
  all.equal(
    colMeans(df_to_test2), 
    colMeans(df_reference2),
    tolerance = 1e-4))
# This should not pass, but since the error I added is centered around 0 it is removed by `colMeans`.

testthat::expect_equal(
  tolerance = 1e-4, # this is relative tolerance
  df_to_test2,
  df_reference2)
#> Error: `df_to_test2` not equal to `df_reference2`.
#> Component "Evap": Mean relative difference: 0.0004877337

^{Created on 2025-02-27 with reprex v2.1.0}

[fabern]

Thanks Mayeul for this major overhaul!
I like that you included regression testing of BiomeE, too.

1. @marcadella Can you confirm that we can safely delete the branches preparation_luluc and preparation_luluc_linked_list ? Former appears to have been merged with a PR here, and latter is indicated to have been merged in a commit here.
2. Note that, I left two code comments regarding the regression testing, suggesting that it should be more strict in the test and more verbose if it fails. (I think they should already show up on this PR thread.)
3. I suggest that this PR should be dealt with first, before we fix the merge conflicts arising from Phydro gw #287

[marcadella]

Up to you. I found it enough to use the mean to detect when I was breaking something. I was also always visualizing some plots.

…

On Thu, Feb 27, 2025, 10:25 Fabian Bernhard ***@***.***> wrote: ***@***.**** commented on this pull request. ------------------------------ In tests/testthat/test-model-runs.R <#282 (comment)>: > + expect_true(all.equal(colMeans(out$data[[1]]$output_daily_tile), colMeans(biomee_p_model_output$data[[1]]$output_daily_tile), tolerance = 1e-4)) + expect_true(all.equal(colMeans(out$data[[1]]$output_annual_tile), colMeans(biomee_p_model_output$data[[1]]$output_annual_tile), tolerance = 1e-4)) + expect_true(all.equal(colMeans(out$data[[1]]$output_annual_cohorts), colMeans(biomee_p_model_output$data[[1]]$output_annual_cohorts), tolerance = 1e-4)) + Illustrating it with the real output, and mocked up errors: library(rsofun) library(dplyr) library(tidyr)df_reference2 <- rsofun::biomee_p_model_output$data[[1]]$output_daily_tile |> tidyr::tibble()df_to_test2 <- df_reference2 |> dplyr::mutate(Evap = Evap * (1 + runif(n(), min = -1e-3, max = 1e-3)), GPP = GPP * (1 + runif(n(), min = -1e-4, max = 1e-4)), slowSOM = slowSOM * (1 + runif(n(), min = -1e-4, max = 1e-4))) testthat::expect_true( all.equal( colMeans(df_to_test2), colMeans(df_reference2), tolerance = 1e-4))# This should not pass, but since the error I added is centered around 0 it is removed by `colMeans`. testthat::expect_equal( tolerance = 1e-4, # this is relative tolerance df_to_test2, df_reference2)#> Error: `df_to_test2` not equal to `df_reference2`.#> Component "Evap": Mean relative difference: 0.0004877337 Created on 2025-02-27 with reprex v2.1.0 <https://reprex.tidyverse.org> — Reply to this email directly, view it on GitHub <#282 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AJMC4KCVCIO3FHVRHNTXQV32R3KX5AVCNFSM6AAAAABWLZ3LWCVHI2DSMVQWIX3LMV43YUDVNRWFEZLROVSXG5CSMV3GSZLXHMZDMNBXGE2DCOBQG4> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[fabern]

The LULUC change appeared to have changed the initial conditions regarding LAI and sapwood to hardwood ratio of carbon. This affects the initial trajectory but appears to smooth out after 8 years.

No action needs to be taken. But it's good to have this documented here.

Detailed observations:
LAI:
The first year output, shows the LAI outputs at only half of the output before the LULUC change.
This links then also to leafC (in line 388).
Sapwood/Hardwood:
Moreover, with the LULUC change, at the end of first year total carbon in sapwood and hardwood is similar, but is distributed differently: sapwood SW_C is only half and hardwood HW_C is double w.r.t the outputs prior to the LULUC change (lines 390 and 391).

These changes are present in the first year ref_BiomeE_Pmodel_odt_yr1 (and might be linked to initialisation or allometry), but disappear for later years in the simulation. Changes to the reference values in ref_BiomeE_Pmodel_odt_yr251 of these quantities are indeed only minor (lines 426 and 427).

Obseravations for ref_BiomeE_Pmodel_oat (and per cohort comparing ref_BiomeE_Pmodel_oac_yr1 vs ref_BiomeE_Pmodel_oac_yr251) are similar. With the LULUC change, in the first two years there is less carbon in SapwoodC than in WoodC, this switches and starting from year 8,9 more carbon in sapwood than in hardwood. By the end of the 250 year simulation the ratio is 4.9/3.7 in favor of sapwood. Prior to the LULUC change this ratio was approximately this from year 1. (lines 469 and 470).

[stineb]

One point @fabern and I discussed:

• Before merging, we should perform a test simulation that includes a spinup and the transient simulation using LUH forcing for a single gridcell, constant CO2, N-deposition, and constant (recycled) climate. We should be looking at the net C balance which should be zero if there was no LUC by definition.