Merge branch 'develop' into phase4-benchmarking-integration

Author: ayushman1210

models/sipnet/NEWS.md

@@ -1,5 +1,16 @@
 # PEcAn.SIPNET 1.10.0.9000
 
+* Improvements to the job.sh written by `write_segmented_configs`:
+  - Now places README.txt, segments.csv, and the full log files from each
+    segment, into the outdir (as was already done for one-segment runs).
+  - Now respects `settings$model$delete.raw` by deleting sipnet.out from segment
+    dirs as well as the from the job outdir
+  - No longer calls model2netcdf.SIPNET() inside each segment dir (it was
+    redundant with the whole-job netcdf output).
+  - Fixed incorrect run paths in the settings$host block passed to single segments
+* `model2netcdf.SIPNET` takes `LAI` from sipnet.out if present (which requires
+    Sipnet > v2.2). If it is not present, LAI is calculated as
+    `plantLeafC / leafCSpWt` as previously.
 * `split_inputs.SIPNET` now avoids internal time format conversions, giving a
   substantial speedup and reduced memory use when processing multi-year files.
 * `model2netcdf.SIPNET` now detects the number of timesteps per day by taking the maximum count across all days in the first simulation year, rather than reading only from day 1. This prevents a factor-of-N error in flux unit conversions when the first day of output is partial (fewer timesteps than a complete day) (#3624, #3989).

models/sipnet/R/model2netcdf.SIPNET.R

@@ -150,16 +150,21 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date,
     )
 
 
-  # calculate LAI for standard output
-  # LAI = plantLeafC / leafCSpWt
-  # both operands are in carbon units (gC/m2 and gC/m2_leaf),
-  # so no carbon fraction conversion (e.g. cFracLeaf) is needed.
-  param <- utils::read.table(file.path(gsub(pattern = "/out/",
-                                            replacement = "/run/", x = outdir),
-                                       "sipnet.param"), stringsAsFactors = FALSE)
-  leafCSpWt <- param[param[, 1] == "leafCSpWt", 2]
-  SLA <- 1000 / leafCSpWt  # m2 leaf / kg C
-
+  if (!("LAI" %in% colnames(sipnet_output))) {
+    # When LAI not reported, calculate it from leaf C and leaf specific weight,
+    # with the latter read from the parameter file.
+    # Note the hardcoded + undocumented assumption that ../run/sipnet.param
+    # exists and contains the values that were used to generate this output.
+    # LAI = plantLeafC [gC/m2] / leafCSpWt [gC/m2 leaf] = m2 leaf / m2
+    # both operands are in carbon units (gC/m2 and gC/m2_leaf),
+    # so no carbon fraction conversion (e.g. cFracLeaf) is needed.
+    param <- utils::read.table(file.path(gsub(pattern = "/out/",
+                                              replacement = "/run/", x = outdir),
+                                         "sipnet.param"), stringsAsFactors = FALSE)
+    leafCSpWt <- param[param[, 1] == "leafCSpWt", 2]
+    SLA <- 1000 / leafCSpWt  # m2 leaf / kg C
+    sipnet_output$LAI = sipnet_output$plantLeafC * SLA
+  }
 
   ### Loop over years in SIPNET output to create separate netCDF outputs
   for (y in year_seq) {
@@ -208,7 +213,7 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date,
       "litter_carbon_content" = sub.sipnet.output$litter,
       "fine_root_carbon_content" = sub.sipnet.output$fineRootC,
       "coarse_root_carbon_content" = sub.sipnet.output$coarseRootC,
-      "LAI" = sub.sipnet.output$plantLeafC * SLA,
+      "LAI" = sub.sipnet.output$LAI,
       "TotLivBiom" = sub.sipnet.output$plantWoodC + sub.sipnet.output$plantLeafC +
                        sub.sipnet.output$coarseRootC + sub.sipnet.output$fineRootC,
       "TotSoilCarb" = sub.sipnet.output$soil + sub.sipnet.output$litter,

models/sipnet/R/write.configs.SIPNET.R

@@ -121,6 +121,16 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
     has_microbeInit = rev_str == "v1"
   )
 
+
+  # find out where to write run/ouput
+  rundir <- file.path(settings$host$rundir, as.character(run.id))
+  outdir <- file.path(settings$host$outdir, as.character(run.id))
+  if (is.null(settings$host$qsub) && (settings$host$name == "localhost")) {
+    rundir <- file.path(settings$rundir, as.character(run.id))
+    outdir <- file.path(settings$modeloutdir, as.character(run.id))
+  }
+
+
   ### WRITE sipnet.in
   template.in <- system.file(
     paste0("sipnet.in_", rev_str),
@@ -140,7 +150,7 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
   }
   config.text <- update_flag_lines(config.text, user_flags)
 
-  writeLines(config.text, con = file.path(settings$rundir, run.id, "sipnet.in"))
+  writeLines(config.text, con = file.path(rundir, "sipnet.in"))
 
   ### WRITE *.clim
   template.clim <- settings$run$inputs$met$path  ## read from settings
@@ -152,14 +162,6 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
   }
   PEcAn.logger::logger.info(paste0("Writing SIPNET configs with input ", template.clim))
 
-  # find out where to write run/ouput
-  rundir <- file.path(settings$host$rundir, as.character(run.id))
-  outdir <- file.path(settings$host$outdir, as.character(run.id))
-  if (is.null(settings$host$qsub) && (settings$host$name == "localhost")) {
-    rundir <- file.path(settings$rundir, as.character(run.id))
-    outdir <- file.path(settings$modeloutdir, as.character(run.id))
-  }
-
   # create launch script (which will create symlink)
   if (!is.null(settings$model$jobtemplate) && file.exists(settings$model$jobtemplate)) {
     jobsh <- readLines(con = settings$model$jobtemplate, n = -1)
@@ -251,8 +253,8 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
   }
   jobsh <- gsub("@DELETE.RAW@", settings$model$delete.raw, jobsh)
 
-  writeLines(jobsh, con = file.path(settings$rundir, run.id, "job.sh"))
-  Sys.chmod(file.path(settings$rundir, run.id, "job.sh"))
+  writeLines(jobsh, con = file.path(rundir, "job.sh"))
+  Sys.chmod(file.path(rundir, "job.sh"))
 
 
   ### Copy event file
@@ -268,7 +270,7 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
   ### WRITE *.param-spatial
   if (caps$has_param_spatial) {
     template.paramSpatial <- system.file("template.param-spatial", package = "PEcAn.SIPNET")
-    file.copy(template.paramSpatial, file.path(settings$rundir, run.id, "sipnet.param-spatial"))
+    file.copy(template.paramSpatial, file.path(rundir, "sipnet.param-spatial"))
   }
 
   ### WRITE *.param
@@ -1028,12 +1030,11 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
     }
 
   }
-  if (file.exists(file.path(settings$rundir, run.id, "sipnet.param"))) {
+  if (file.exists(file.path(rundir, "sipnet.param"))) {
     file.rename(
-      file.path(settings$rundir, run.id, "sipnet.param"),
+      file.path(rundir, "sipnet.param"),
       file.path(
-        settings$rundir,
-        run.id,
+        rundir,
         paste0("sipnet_", lubridate::year(settings$run$start.date), "_", lubridate::year(settings$run$end.date), ".param")
       )
     )
@@ -1042,7 +1043,7 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
 
   utils::write.table(
     param,
-    file.path(settings$rundir, run.id, "sipnet.param"),
+    file.path(rundir, "sipnet.param"),
     row.names = FALSE,
     col.names = FALSE,
     quote = FALSE

models/sipnet/R/write_segmented_configs.R

@@ -97,6 +97,15 @@ write_segmented_configs.SIPNET <- function(settings, input_design = NULL, ...) {
         by = "ens_num",
         relationship = "many-to-one")
   }
+  # Use a custom job.sh template with no model2netcdf step;
+  # we convert all segments together at the end of the timeseries instead.
+  # ...unless user provided their own template, in which case use it as-is
+  if (is.null(settings$model$jobtemplate) || !file.exists(settings$model$jobtemplate)) {
+    settings$model$jobtemplate <- system.file(
+      "template_singlesegment.job",
+      package = "PEcAn.SIPNET"
+    )
+  }
 
   new_jobfiles <- character()
 
@@ -218,6 +227,7 @@ write_segment_configs <- function(
   )
 
   jobsh_files <- character()
+  job_logfiles <- character()
 
   for (isegment in seq_len(nrow(segments))) {
     segment <- segments[isegment, ]
@@ -259,6 +269,18 @@ write_segment_configs <- function(
     segment_settings[["outdir"]] <- segment_outdir
     segment_settings[["modeloutdir"]] <- segment_outdir
     segment_settings[["rundir"]] <- segment_rundir
+    segment_settings[[c("host", "rundir")]] <-  sub(
+      pattern = settings$rundir,
+      replacement = settings$host$rundir,
+      x = segment_rundir,
+      fixed = TRUE
+    )
+    segment_settings[[c("host", "outdir")]] <-  sub(
+      pattern = settings$modeloutdir,
+      replacement = settings$host$outdir,
+      x = segment_outdir,
+      fixed = TRUE
+    )
     segment_settings[[c("run", "start.date")]] <- dstart
     segment_settings[[c("run", "end.date")]] <- dend
     segment_settings[[c("run", "inputs")]] <- segment_inputs
@@ -291,12 +313,19 @@ write_segment_configs <- function(
       run.id = runid_dummy
     )
 
+    segment_log <- file.path(segment_settings$modeloutdir,
+                             runid_dummy,
+                             "logfile.txt")
     segment_jobsh <- file.path(segment_settings$rundir, runid_dummy, "job.sh")
     stopifnot(file.exists(segment_jobsh))
+    job_logfiles <- c(job_logfiles, segment_log)
     jobsh_files <- c(jobsh_files, segment_jobsh)
   }
 
   # Now, get the run's jobsh file
+  # NB this doesn't do any path expansion or @TEMPLATE@ string replacement --
+  # those are done in the per-segment job.sh.
+  # Note especially that means no @HOST_SETUP@ / @HOST_TEARDOWN@ in this script
   run_jobsh <- file.path(run_dir, "job.sh")
   target_sipnet_out <- file.path(run_modeloutdir, "sipnet.out")
   segmented_jobsh_file <- file.path(run_dir, "job_segmented.sh")
@@ -326,10 +355,31 @@ write_segment_configs <- function(
         sprintf("sitelon = %s", as.character(settings$run$site$lon)),
         sprintf("start_date = %s", shQuote(settings$run$start.date)),
         sprintf("end_date = %s", shQuote(settings$run$end.date)),
+        sprintf("delete.raw = %s", shQuote(settings$model$delete.raw)),
         sprintf("revision = %s", shQuote(settings$model$revision)),
         sep = ", "
       )
-    )
+    ),
+    "",
+    "# copy readme with specs to output",
+    paste("cp", file.path(run_dir, "README.txt"), file.path(run_modeloutdir, "README.txt")),
+    paste("cp", file.path(run_dir, "segments.csv"), file.path(run_modeloutdir, "segments.csv")),
+    "",
+    "# Concatenate segment log files",
+    paste("echo \"\n--> contents of", job_logfiles, ":\" && cat", job_logfiles),
+    "",
+    if (isTRUE(as.logical(settings$model$delete.raw))) {
+      # Removing all outputs that have been safely copied to the outdir
+      c(
+        "# Remove per-segment outputs & logs after concatenating to job outdir",
+        sprintf("find %s -name sipnet.out -delete", segment_rootdir),
+        sprintf("find %s -name logfile.txt -delete", segment_rootdir),
+        sprintf("rm %s/README.txt", run_dir),
+        sprintf("rm %s/segments.csv", run_dir)
+      )
+    },
+    "",
+    "echo -e \"MODEL FINISHED\nLogfile is located at '${OUTDIR}/logfile.txt'\" >&3"
   )
   writeLines(segmented_jobsh_lines, segmented_jobsh_file)
   if (replace_and_link) {

models/sipnet/inst/template_singlesegment.job

@@ -0,0 +1,53 @@
+#!/bin/bash
+
+# redirect output
+exec 3>&1
+exec &> "@OUTDIR@/logfile.txt"
+
+# host specific setup
+@HOST_SETUP@
+
+# cdo setup
+@CDO_SETUP@
+
+# create output folder
+mkdir -p "@OUTDIR@"
+
+# Convert any relative paths to absolute
+# (otherwise we'd lose track of them when cd'ing into rundir)
+OUTDIR=$(cd "@OUTDIR@" && pwd -P)
+RUNDIR=$(cd "@RUNDIR@" && pwd -P)
+SITE_MET=$(cd $(dirname "@SITE_MET@") && pwd -P)/$(basename "@SITE_MET@")
+BINARY=$(cd $(dirname "@BINARY@") && pwd -P)/$(basename "@BINARY@")
+
+# see if application needs running
+if [ ! -e "${OUTDIR}/sipnet.out" ]; then
+  cd "$RUNDIR"
+  ln -s "${SITE_MET}" sipnet.clim
+
+  "${BINARY}"
+  STATUS=$?
+  
+  # copy output
+  mv "${RUNDIR}/sipnet.out" "$OUTDIR"
+
+  # check the status
+  if [ $STATUS -ne 0 ]; then
+    echo -e "ERROR IN MODEL RUN\nLogfile is located at '${OUTDIR}/logfile.txt'" >&3
+    exit $STATUS
+  fi
+
+  # conversion to MsTMIP not done in this single-segment step.
+  # The whole-job script will handle conversion after all segments have run.
+fi
+
+# copy readme with specs to output
+cp  "${RUNDIR}/README.txt" "${OUTDIR}/README.txt"
+
+# run getdata to extract right variables
+
+# host specific teardown
+@HOST_TEARDOWN@
+
+# all done
+echo -e "MODEL FINISHED\nLogfile is located at '${OUTDIR}/logfile.txt'" >&3

models/sipnet/tests/testthat/test-model2netcdf.SIPNET.R

@@ -291,3 +291,22 @@ test_that("out_day is correct when the first day of output is partial", {
   ts_correct <- 86400 / 2
   expect_equal(gpp, rep(base_row$gpp * 1e-3 / ts_correct, 5), tolerance = 1e-10)
 })
+
+test_that("LAI taken from leafCSpWt if not present in output", {
+  dat <- make_v2_sipnet(2)
+  # setup_sipnet_test sets leafCSpWt = 32 => expect LAI of 1 and 10
+  dat$plantLeafC = c(32, 320)
+
+  # First without LAI in dat => looks for sipnet.param in run dir
+  paths <- setup_sipnet_test(dat, notes_line = NULL)
+  nc1 <- ncdf4::nc_open(file.path(paths$outdir, "2002.nc"))
+  on.exit(ncdf4::nc_close(nc1), add = TRUE)
+  expect_equal(as.vector(ncdf4::ncvar_get(nc1, "LAI")), c(1, 10))
+
+  # Now with LAI specified => uses it as-is
+  dat$LAI = c(2, 3)
+  paths <- setup_sipnet_test(dat, notes_line = NULL)
+  nc2 <- ncdf4::nc_open(file.path(paths$outdir, "2002.nc"))
+  on.exit(ncdf4::nc_close(nc2), add = TRUE)
+  expect_equal(as.vector(ncdf4::ncvar_get(nc2, "LAI")), c(2, 3))
+})