SIP400 Address harvest and leaf off issues

docs/CHANGELOG.md

@@ -26,6 +26,8 @@ sections to include in release notes:
 ### Added
 
 ### Fixed
+- Terminating harvest events will now remove all biomass at end of time step, instead of allowing small amounts to remain in the plant pools due to growth in that step. (#400)
+- Leaf-off mechanics refined to ensure that leaf-off events do not remove more leaf carbon than is available in the leaf pool. (#400)
 
 ### Changed
 - SIPNET will now error instead of warning when an environment pool goes negative. This is a change from previous behavior where SIPNET would log a warning and continue running. (#398)

src/common/exitCodes.h

@@ -23,7 +23,8 @@ typedef enum {
   EXIT_CODE_FILE_OPEN_OR_READ_ERROR = 6,
   EXIT_CODE_INTERNAL_ERROR = 7,
   EXIT_CODE_BAD_CLI_ARGUMENT = 8,
-  EXIT_CODE_BAD_RESTART_PARAMETER = 9
+  EXIT_CODE_BAD_RESTART_PARAMETER = 9,
+  EXIT_CODE_MEMORY_ALLOCATION_FAILURE = 10
 } exit_code_t;
 
 #endif

src/common/util.c

@@ -76,3 +76,109 @@ double calcRatio(const double num, const double den) {
 
 // For global linkage
 extern inline double unitClip(double preClip);
+
+DynamicString *dsCreate(size_t initial_capacity) {
+  DynamicString *ds = malloc(sizeof(DynamicString));
+  if (!ds)
+    return NULL;
+
+  // Ensure we have room for at least a null terminator
+  ds->capacity = (initial_capacity > 0) ? initial_capacity : 16;
+  ds->buffer = malloc(ds->capacity * sizeof(char));
+
+  if (!ds->buffer) {
+    free(ds);
+    return NULL;
+  }
+
+  ds->buffer[0] = '\0';  // Start with an empty string
+  ds->length = 0;
+  return ds;
+}
+
+int dsAppend(DynamicString *ds, const char *str) {
+  if (!ds || !str)
+    return 0;
+
+  size_t append_len = strlen(str);
+  // +1 is crucial to ensure room for the null terminator
+  size_t needed_capacity = ds->length + append_len + 1;
+
+  // Double capacity until it's large enough for the new content
+  if (needed_capacity > ds->capacity) {
+    size_t new_capacity = ds->capacity * 2;
+    while (new_capacity < needed_capacity) {
+      new_capacity *= 2;
+    }
+
+    // Safely reallocate using a temporary pointer
+    char *temp = realloc(ds->buffer, new_capacity * sizeof(char));
+    if (!temp) {
+      return 0;  // Reallocation failed, original data remains intact
+    }
+
+    ds->buffer = temp;
+    ds->capacity = new_capacity;
+  }
+
+  // Copy the new string over the old null terminator
+  memcpy(ds->buffer + ds->length, str, append_len);
+  ds->length += append_len;
+  ds->buffer[ds->length] = '\0';  // Manually place the new null terminator
+
+  return 1;  // Success
+}
+
+int dsAppendFormatted(DynamicString *ds, const char *format, ...) {
+  if (!ds || !format)
+    return 0;
+
+  // 1. Determine how much space the formatted text needs
+  va_list args;
+  va_start(args, format);
+  // Make a copy of args because vsnprintf consumes the list
+  va_list args_copy;
+  va_copy(args_copy, args);
+
+  // Pass NULL and 0 to just count the characters needed
+  int formatted_len = vsnprintf(NULL, 0, format, args_copy);
+  va_end(args_copy);
+
+  if (formatted_len < 0) {
+    va_end(args);
+    return 0;  // Formatting error occurred
+  }
+
+  // 2. Ensure the buffer is big enough
+  size_t needed_capacity = ds->length + (size_t)formatted_len + 1;
+  if (needed_capacity > ds->capacity) {
+    size_t new_capacity = ds->capacity * 2;
+    while (new_capacity < needed_capacity) {
+      new_capacity *= 2;
+    }
+
+    char *temp = realloc(ds->buffer, new_capacity * sizeof(char));
+    if (!temp) {
+      va_end(args);
+      return 0;  // Reallocation failed
+    }
+    ds->buffer = temp;
+    ds->capacity = new_capacity;
+  }
+
+  // 3. Write the formatted string directly into the builder's buffer
+  // Write directly to the position of the current null terminator
+  vsnprintf(ds->buffer + ds->length, (size_t)formatted_len + 1, format, args);
+  va_end(args);
+
+  // 4. Update the length of the string builder
+  ds->length += (size_t)formatted_len;
+  return 1;
+}
+
+void dsFree(DynamicString *ds) {
+  if (ds) {
+    free(ds->buffer);
+    free(ds);
+  }
+}

src/common/util.h

@@ -9,7 +9,10 @@
 #define UTIL_H
 
 #include <math.h>
+#include <stdarg.h>  // Required for va_list, va_start, va_end
 #include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
 
 #define TINY 0.000001  // to avoid those nasty divide-by-zero errors
 
@@ -37,4 +40,23 @@ double calcRatio(double num, double den);
  */
 inline double unitClip(double preClip) { return fmin(fmax(preClip, 0.0), 1.0); }
 
+// DYNAMIC STRING
+typedef struct DynamicStringStruct {
+  char *buffer;  // Pointer to the character array
+  size_t length;  // Number of characters currently in the string
+  size_t capacity;  // Total allocated space (including null terminator)
+} DynamicString;
+
+// Initialize the builder with a reasonable default capacity
+DynamicString *dsCreate(size_t initial_capacity);
+
+// Append text, automatically doubling capacity if needed
+int dsAppend(DynamicString *ds, const char *str);
+
+// Append formatted text using printf-style syntax
+int dsAppendFormatted(DynamicString *ds, const char *format, ...);
+
+// Free all memory associated with the builder
+void dsFree(DynamicString *ds);
+
 #endif

src/sipnet/balance.c

@@ -47,7 +47,8 @@ void updateBalanceTrackerPostClamp(void) {
   balanceTracker.clampedC = balanceTracker.finalC - balanceTracker.postTotalC;
   if (balanceTracker.clampedC < -EPS) {
     // This shouldn't happen, by construction
-    logInternalError("Non-negative clamping has cause carbon loss\n");
+    logInternalError("Non-negative clamping has cause carbon loss %f\n",
+                     balanceTracker.clampedC);
   }
   if (balanceTracker.clampedC < EPS) {
     balanceTracker.clampedC = 0;
@@ -56,7 +57,8 @@ void updateBalanceTrackerPostClamp(void) {
   balanceTracker.clampedN = balanceTracker.finalN - balanceTracker.postTotalN;
   if (balanceTracker.clampedN < -EPS) {
     // This shouldn't happen, by construction
-    logInternalError("Non-negative clamping has cause nitrogen loss\n");
+    logInternalError("Non-negative clamping has cause nitrogen loss %f\n",
+                     balanceTracker.clampedN);
   }
   if (balanceTracker.clampedN < EPS) {
     balanceTracker.clampedN = 0;

src/sipnet/debug_log.c

@@ -4,11 +4,44 @@
 
 #include "debug_log.h"
 
+#include "events.h"
 #include "common/exitCodes.h"
 #include "common/logging.h"
 #include "common/context.h"
 #include "common/util.h"
 #include "state.h"
+
+#define NUM_LOGGED_ENVI_FIELDS 13
+#define NUM_LOGGED_FLUX_FIELDS 59
+#define NUM_LOGGED_TRACKER_FIELDS 33
+#define NUM_LOGGED_PHEN_TRACKER_FIELDS 3
+#define NUM_LOGGED_SURVIVAL_FIELDS 1
+#define NUM_LOGGED_EVENT_TRACKER_FIELDS 7
+
+#define DEBUG_LAYOUT_ENVI_SIZE (8 * NUM_LOGGED_ENVI_FIELDS)
+#define DEBUG_LAYOUT_FLUX_SIZE (8 * NUM_LOGGED_FLUX_FIELDS)
+// The Trackers struct is not all doubles, but the int(s) get padded to 8 bytes
+#define DEBUG_LAYOUT_TRACKER_SIZE (8 * NUM_LOGGED_TRACKER_FIELDS)
+#define DEBUG_LAYOUT_PHEN_SIZE (4 * NUM_LOGGED_PHEN_TRACKER_FIELDS)
+#define DEBUG_LAYOUT_SURVIVAL_SIZE (4 * NUM_LOGGED_SURVIVAL_FIELDS)
+#define DEBUG_LAYOUT_EVENT_SIZE (8 * NUM_LOGGED_EVENT_TRACKER_FIELDS)
+
+_Static_assert(sizeof(Envi) == DEBUG_LAYOUT_ENVI_SIZE,
+               "Debug log schema drift: Envi changed; update debug_log.c");
+_Static_assert(sizeof(Fluxes) == DEBUG_LAYOUT_FLUX_SIZE,
+               "Debug log schema drift: Fluxes changed; update debug_log.c");
+_Static_assert(sizeof(Trackers) == DEBUG_LAYOUT_TRACKER_SIZE,
+               "Debug log schema drift: Trackers changed; update debug_log.c");
+_Static_assert(
+    sizeof(PhenologyTrackers) == DEBUG_LAYOUT_PHEN_SIZE,
+    "Debug log schema drift: PhenologyTrackers changed; update debug_log.c");
+_Static_assert(
+    sizeof(PlantSurvivalTracker) == DEBUG_LAYOUT_SURVIVAL_SIZE,
+    "Debug log schema drift: PlantSurvival changed; update debug_log.c");
+_Static_assert(
+    sizeof(EventTrackers) == DEBUG_LAYOUT_EVENT_SIZE,
+    "Debug log schema drift: EventTrackers changed; update debug_log.c");
+
 typedef enum DebugFieldType {
   DEBUG_FIELD_INT = 0,
   DEBUG_FIELD_DOUBLE = 1
@@ -20,18 +53,13 @@ typedef struct DebugField {
   const void *value;
 } DebugField;
 
-#define NUM_LOGGED_ENVI_FIELDS 13
-#define NUM_LOGGED_FLUX_FIELDS 56
-#define NUM_LOGGED_TRACKER_FIELDS 33
-#define NUM_LOGGED_PHEN_TRACKER_FIELDS 3
-#define NUM_LOGGED_SURVIVAL_FIELDS 1
-
 typedef struct DebugFieldArrays {
   DebugField enviDF[NUM_LOGGED_ENVI_FIELDS];
   DebugField fluxDF[NUM_LOGGED_FLUX_FIELDS];
   DebugField trackerDF[NUM_LOGGED_TRACKER_FIELDS];
   DebugField phenoDF[NUM_LOGGED_PHEN_TRACKER_FIELDS];
   DebugField survivalDF[NUM_LOGGED_SURVIVAL_FIELDS];
+  DebugField eventDF[NUM_LOGGED_EVENT_TRACKER_FIELDS];
 } DebugFieldArrays;
 
 static DebugFieldArrays *debugFields = NULL;
@@ -61,7 +89,11 @@ void initDebugArrays() {
   debugFields->enviDF[ind++] = (DebugField){"soilOrgN", DEBUG_FIELD_DOUBLE, &envi.soilOrgN},
   debugFields->enviDF[ind++] = (DebugField){"litterN", DEBUG_FIELD_DOUBLE, &envi.litterN},
   debugFields->enviDF[ind++] = (DebugField){"plantStorageN", DEBUG_FIELD_DOUBLE, &envi.plantStorageN},
-  debugFields->enviDF[ind  ] = (DebugField){"plantCAccountingDelta", DEBUG_FIELD_DOUBLE,&envi.plantCAccountingDelta};
+  debugFields->enviDF[ind++] = (DebugField){"plantCAccountingDelta", DEBUG_FIELD_DOUBLE,&envi.plantCAccountingDelta};
+  if (ind != NUM_LOGGED_ENVI_FIELDS) {
+    logInternalError("Debug log array size mismatch: enviDF\n");
+    exit(EXIT_CODE_INTERNAL_ERROR);
+  }
 
   ind = 0;
   debugFields->fluxDF[ind++] = (DebugField){"photosynthesis", DEBUG_FIELD_DOUBLE, &fluxes.photosynthesis};
@@ -90,6 +122,7 @@ void initDebugArrays() {
   debugFields->fluxDF[ind++] = (DebugField){"woodCreation", DEBUG_FIELD_DOUBLE, &fluxes.woodCreation},
   debugFields->fluxDF[ind++] = (DebugField){"leafOnCreation", DEBUG_FIELD_DOUBLE, &fluxes.leafOnCreation},
   debugFields->fluxDF[ind++] = (DebugField){"leafOnCreationFromWood", DEBUG_FIELD_DOUBLE, &fluxes.leafOnCreationFromWood},
+  debugFields->fluxDF[ind++] = (DebugField){"leafOffLitter", DEBUG_FIELD_DOUBLE, &fluxes.leafOffLitter},
   debugFields->fluxDF[ind++] = (DebugField){"nVolatilization", DEBUG_FIELD_DOUBLE, &fluxes.nVolatilization},
   debugFields->fluxDF[ind++] = (DebugField){"nLeaching", DEBUG_FIELD_DOUBLE, &fluxes.nLeaching},
   debugFields->fluxDF[ind++] = (DebugField){"nOrgSoil", DEBUG_FIELD_DOUBLE, &fluxes.nOrgSoil},
@@ -101,6 +134,7 @@ void initDebugArrays() {
   debugFields->fluxDF[ind++] = (DebugField){"reductionNResorption", DEBUG_FIELD_DOUBLE, &fluxes.reductionNResorption},
   debugFields->fluxDF[ind++] = (DebugField){"eventLeafC", DEBUG_FIELD_DOUBLE, &fluxes.eventLeafC},
   debugFields->fluxDF[ind++] = (DebugField){"eventWoodC", DEBUG_FIELD_DOUBLE, &fluxes.eventWoodC},
+  debugFields->fluxDF[ind++] = (DebugField){"eventAccountingC", DEBUG_FIELD_DOUBLE, &fluxes.eventAccountingC},
   debugFields->fluxDF[ind++] = (DebugField){"eventFineRootC", DEBUG_FIELD_DOUBLE, &fluxes.eventFineRootC},
   debugFields->fluxDF[ind++] = (DebugField){"eventCoarseRootC", DEBUG_FIELD_DOUBLE, &fluxes.eventCoarseRootC},
   debugFields->fluxDF[ind++] = (DebugField){"eventEvap", DEBUG_FIELD_DOUBLE, &fluxes.eventEvap},
@@ -116,10 +150,15 @@ void initDebugArrays() {
   debugFields->fluxDF[ind++] = (DebugField){"eventOutputN", DEBUG_FIELD_DOUBLE, &fluxes.eventOutputN},
   debugFields->fluxDF[ind++] = (DebugField){"eventLeafOnCreation", DEBUG_FIELD_DOUBLE, &fluxes.eventLeafOnCreation},
   debugFields->fluxDF[ind++] = (DebugField){"eventLeafOnCreationFromWood", DEBUG_FIELD_DOUBLE, &fluxes.eventLeafOnCreationFromWood},
-  debugFields->fluxDF[ind++] = (DebugField){"eventLeafOffLitter", DEBUG_FIELD_DOUBLE, &fluxes.eventLeafOffLitter},
+  debugFields->fluxDF[ind++] = (DebugField){"eventLeafOffLitterC", DEBUG_FIELD_DOUBLE, &fluxes.eventLeafOffLitterC},
+  debugFields->fluxDF[ind++] = (DebugField){"eventLeafOffLitterN", DEBUG_FIELD_DOUBLE, &fluxes.eventLeafOffLitterN},
   debugFields->fluxDF[ind++] = (DebugField){"eventLeafOffNResorption", DEBUG_FIELD_DOUBLE, &fluxes.eventLeafOffNResorption},
   debugFields->fluxDF[ind++] = (DebugField){"soilMethane", DEBUG_FIELD_DOUBLE, &fluxes.soilMethane},
-  debugFields->fluxDF[ind  ] = (DebugField){"litterMethane", DEBUG_FIELD_DOUBLE, &fluxes.litterMethane};
+  debugFields->fluxDF[ind++] = (DebugField){"litterMethane", DEBUG_FIELD_DOUBLE, &fluxes.litterMethane};
+  if (ind != NUM_LOGGED_FLUX_FIELDS) {
+    logInternalError("Debug log array size mismatch: fluxDF\n");
+    exit(EXIT_CODE_INTERNAL_ERROR);
+  }
 
   ind = 0;
   debugFields->trackerDF[ind++] = (DebugField){"gpp", DEBUG_FIELD_DOUBLE, &trackers.gpp},
@@ -154,15 +193,40 @@ void initDebugArrays() {
   debugFields->trackerDF[ind++] = (DebugField){"nLeaching", DEBUG_FIELD_DOUBLE, &trackers.nLeaching},
   debugFields->trackerDF[ind++] = (DebugField){"nFixation", DEBUG_FIELD_DOUBLE, &trackers.nFixation},
   debugFields->trackerDF[ind++] = (DebugField){"nUptake", DEBUG_FIELD_DOUBLE, &trackers.nUptake};
-  debugFields->trackerDF[ind  ] = (DebugField){"meanNPP", DEBUG_FIELD_DOUBLE, &trackers.meanNPP};
+  debugFields->trackerDF[ind++] = (DebugField){"meanNPP", DEBUG_FIELD_DOUBLE, &trackers.meanNPP};
+  if (ind != NUM_LOGGED_TRACKER_FIELDS) {
+    logInternalError("Debug log array size mismatch: trackerDF\n");
+    exit(EXIT_CODE_INTERNAL_ERROR);
+  }
 
   ind = 0;
   debugFields->phenoDF[ind++] = (DebugField){"didLeafGrowth", DEBUG_FIELD_INT, &phenologyTrackers.didLeafGrowth},
   debugFields->phenoDF[ind++] = (DebugField){"didLeafFall", DEBUG_FIELD_INT, &phenologyTrackers.didLeafFall},
-  debugFields->phenoDF[ind  ] = (DebugField){"lastYear", DEBUG_FIELD_INT, &phenologyTrackers.lastYear};
+  debugFields->phenoDF[ind++] = (DebugField){"lastYear", DEBUG_FIELD_INT, &phenologyTrackers.lastYear};
+  if (ind != NUM_LOGGED_PHEN_TRACKER_FIELDS) {
+    logInternalError("Debug log array size mismatch: phenoDF\n");
+    exit(EXIT_CODE_INTERNAL_ERROR);
+  }
 
   ind = 0;
-  debugFields->survivalDF[ind] = (DebugField){"isAlive", DEBUG_FIELD_INT, &plantSurvivalTracker.isAlive};
+  debugFields->survivalDF[ind++] = (DebugField){"isAlive", DEBUG_FIELD_INT, &plantSurvivalTracker.isAlive};
+  if (ind != NUM_LOGGED_SURVIVAL_FIELDS) {
+    logInternalError("Debug log array size mismatch: survivalDF\n");
+    exit(EXIT_CODE_INTERNAL_ERROR);
+  }
+
+  ind = 0;
+  debugFields->eventDF[ind++] = (DebugField){"d_till_mod", DEBUG_FIELD_DOUBLE, &eventTrackers.d_till_mod};
+  debugFields->eventDF[ind++] = (DebugField){"ht.totalFracRemoved", DEBUG_FIELD_DOUBLE, &eventTrackers.harvestTrackers.totalFracRemoved};
+  debugFields->eventDF[ind++] = (DebugField){"ht.totalFracTransferred", DEBUG_FIELD_DOUBLE, &eventTrackers.harvestTrackers.totalFracTransferred};
+  debugFields->eventDF[ind++] = (DebugField){"ht.totalFracRemovedAbove", DEBUG_FIELD_DOUBLE, &eventTrackers.harvestTrackers.totalFracRemovedAbove};
+  debugFields->eventDF[ind++] = (DebugField){"ht.totalFracRemovedBelow", DEBUG_FIELD_DOUBLE, &eventTrackers.harvestTrackers.totalFracRemovedBelow};
+  debugFields->eventDF[ind++] = (DebugField){"ht.totalFracTransferredAbove", DEBUG_FIELD_DOUBLE, &eventTrackers.harvestTrackers.totalFracTransferredAbove};
+  debugFields->eventDF[ind++] = (DebugField){"ht.totalFracTransferredBelow", DEBUG_FIELD_DOUBLE, &eventTrackers.harvestTrackers.totalFracTransferredBelow};
+  if (ind != NUM_LOGGED_EVENT_TRACKER_FIELDS) {
+    logInternalError("Debug log array size mismatch: eventDF\n");
+    exit(EXIT_CODE_INTERNAL_ERROR);
+  }
   // clang-format on
 }
 
@@ -278,6 +342,8 @@ void outputDebugHeaders(DebugLogFiles *debugLogFiles) {
     outputDebugFieldHeader(debugLogFiles->trackers, "s.",
                            debugFields->survivalDF, NUM_LOGGED_SURVIVAL_FIELDS,
                            0);
+    outputDebugFieldHeader(debugLogFiles->trackers, "et.", debugFields->eventDF,
+                           NUM_LOGGED_EVENT_TRACKER_FIELDS, 0);
     fprintf(debugLogFiles->trackers, "\n");
   }
 }
@@ -307,6 +373,9 @@ void outputDebugState(DebugLogFiles *debugLogFiles, int year, int day,
     outputDebugFieldValues(debugLogFiles->trackers, year, day, time,
                            debugFields->survivalDF, NUM_LOGGED_SURVIVAL_FIELDS,
                            0);
+    outputDebugFieldValues(debugLogFiles->trackers, year, day, time,
+                           debugFields->eventDF,
+                           NUM_LOGGED_EVENT_TRACKER_FIELDS, 0);
     fprintf(debugLogFiles->trackers, "\n");
   }
 }