[#416] Remove clamping from balance checks

docs/CHANGELOG.md

@@ -43,6 +43,7 @@ sections to include in release notes:
 ### 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)
 - SIPNET now errors and exits when a mass balance check failure occurs, instead of logging a warning and continuing (#414)
+- Removed the stock clamping adjustment from the mass balance check, which now allows a tolerance of `10 * EPS` (#425)
 
 ### Removed
 

src/sipnet/balance.c

@@ -37,32 +37,6 @@ void updateBalanceTrackerPreUpdate(void) {
 void updateBalanceTrackerPostUpdate(void) {
   // Set the post-update pool totals
   getMassTotals(&balanceTracker.postTotalC, &balanceTracker.postTotalN);
-}
-
-void updateBalanceTrackerPostClamp(void) {
-  getMassTotals(&balanceTracker.finalC, &balanceTracker.finalN);
-
-  // Difference between post-update and post-clamp is the amount we "gained" by
-  // setting negative values to zero
-  balanceTracker.clampedC = balanceTracker.finalC - balanceTracker.postTotalC;
-  if (balanceTracker.clampedC < -EPS) {
-    // This shouldn't happen, by construction
-    logInternalError("Non-negative clamping has caused carbon loss %f\n",
-                     balanceTracker.clampedC);
-  }
-  if (balanceTracker.clampedC < EPS) {
-    balanceTracker.clampedC = 0;
-  }
-
-  balanceTracker.clampedN = balanceTracker.finalN - balanceTracker.postTotalN;
-  if (balanceTracker.clampedN < -EPS) {
-    // This shouldn't happen, by construction
-    logInternalError("Non-negative clamping has caused nitrogen loss %f\n",
-                     balanceTracker.clampedN);
-  }
-  if (balanceTracker.clampedN < EPS) {
-    balanceTracker.clampedN = 0;
-  }
 
   // Calculate the system inputs and outputs
   // CARBON
@@ -90,25 +64,15 @@ void updateBalanceTrackerPostClamp(void) {
     balanceTracker.inputsN *= climate->length;
     balanceTracker.outputsN *= climate->length;
   }
-
-  // Account for gains from clamping
-  balanceTracker.inputsC += balanceTracker.clampedC;
-  if (ctx.nitrogenCycle) {
-    balanceTracker.inputsN += balanceTracker.clampedN;
-  }
 }
 
 void initBalanceTracker(void) {
   // Initialize all to zero
   // Pools
   balanceTracker.preTotalC = 0.0;
   balanceTracker.postTotalC = 0.0;
-  balanceTracker.clampedC = 0.0;
-  balanceTracker.finalC = 0.0;
   balanceTracker.preTotalN = 0.0;
   balanceTracker.postTotalN = 0.0;
-  balanceTracker.clampedN = 0.0;
-  balanceTracker.finalN = 0.0;
 
   // System I/O
   balanceTracker.inputsC = 0.0;
@@ -122,65 +86,60 @@ void initBalanceTracker(void) {
 }
 
 int balanceResidualFailed(double residual) {
-  return !isfinite(residual) || fabs(residual) >= EPS;
+  return !isfinite(residual) || fabs(residual) >= BALANCE_TOLERANCE;
 }
 
 // %.17g round-trips a double, so differences near EPS remain visible even
 // when pool totals are large
 static void logBalanceFailure(const char *element, double residual,
-                              double preTotal, double postTotal,
-                              double finalTotal, double inputs, double outputs,
-                              double clamped) {
+                              double preTotal, double postTotal, double inputs,
+                              double outputs) {
   logInternalError("%s balance check failed: residual=%.17g (tolerance %g) at "
                    "year %d day %d time %.17g\n",
-                   element, residual, EPS, climate->year, climate->day,
-                   climate->time);
-  logInternalError("  %s pools: preTotal=%.17g postUpdateTotal=%.17g "
-                   "finalTotal=%.17g\n",
-                   element, preTotal, postTotal, finalTotal);
-  logInternalError("  %s system: inputs=%.17g (includes clamped) outputs=%.17g "
-                   "clamped=%.17g\n",
-                   element, inputs, outputs, clamped);
+                   element, residual, BALANCE_TOLERANCE, climate->year,
+                   climate->day, climate->time);
+  logInternalError("  %s pools: preTotal=%.17g postTotal=%.17g\n", element,
+                   preTotal, postTotal);
+  logInternalError("  %s system: inputs=%.17g outputs=%.17g\n", element, inputs,
+                   outputs);
 }
 
 void checkBalance(void) {
   // CARBON
   // Pool delta
-  double poolCDelta = balanceTracker.finalC - balanceTracker.preTotalC;
+  double poolCDelta = balanceTracker.postTotalC - balanceTracker.preTotalC;
   // System delta
   double systemCDelta = balanceTracker.inputsC - balanceTracker.outputsC;
   balanceTracker.deltaC = poolCDelta - systemCDelta;
 
   // NITROGEN
   // Pool delta
-  double poolNDelta = balanceTracker.finalN - balanceTracker.preTotalN;
+  double poolNDelta = balanceTracker.postTotalN - balanceTracker.preTotalN;
   // System delta
   double systemNDelta = balanceTracker.outputsN - balanceTracker.inputsN;
   balanceTracker.deltaN = poolNDelta + systemNDelta;
 
-  // Finite residuals smaller than EPS are numerical noise; zero them to avoid
+  // Finite residuals within tolerance are numerical noise; zero them to avoid
   // weird negative-zero issues. NaN and infinity fail this test and are kept.
-  if (fabs(balanceTracker.deltaC) < EPS) {
+  if (fabs(balanceTracker.deltaC) < BALANCE_TOLERANCE) {
     balanceTracker.deltaC = 0.0;
   }
-  if (fabs(balanceTracker.deltaN) < EPS) {
+  if (fabs(balanceTracker.deltaN) < BALANCE_TOLERANCE) {
     balanceTracker.deltaN = 0.0;
   }
 
   int err = 0;
   if (balanceResidualFailed(balanceTracker.deltaC)) {
     err = 1;
     logBalanceFailure("Carbon", balanceTracker.deltaC, balanceTracker.preTotalC,
-                      balanceTracker.postTotalC, balanceTracker.finalC,
-                      balanceTracker.inputsC, balanceTracker.outputsC,
-                      balanceTracker.clampedC);
+                      balanceTracker.postTotalC, balanceTracker.inputsC,
+                      balanceTracker.outputsC);
   }
   if (balanceResidualFailed(balanceTracker.deltaN)) {
     err = 1;
     logBalanceFailure("Nitrogen", balanceTracker.deltaN,
                       balanceTracker.preTotalN, balanceTracker.postTotalN,
-                      balanceTracker.finalN, balanceTracker.inputsN,
-                      balanceTracker.outputsN, balanceTracker.clampedN);
+                      balanceTracker.inputsN, balanceTracker.outputsN);
   }
   if (err) {
     logInternalError("Exiting\n");

src/sipnet/balance.h

@@ -1,10 +1,15 @@
 #ifndef BALANCE_H
 #define BALANCE_H
 
-// Floating-point precision threshold for balance checks and clamping warnings.
-// Differences smaller than this are treated as numerical noise.
+// Floating-point precision threshold. Differences smaller than this are
+// treated as numerical noise.
 #define EPS 1e-8
 
+// Tolerance for the mass balance residual. ensureNonNegativeStocks() may zero
+// several stocks in [-EPS, 0) in a single step, and those small adjustments all
+// end up in the balance residual, so allow for some accumulation.
+#define BALANCE_TOLERANCE (10 * EPS)
+
 typedef struct BalanceTrackerStruct {
   // Mass balance checks:
   //   X_t - X_(t-1) = inputs - outputs + tolerance
@@ -17,23 +22,15 @@ typedef struct BalanceTrackerStruct {
   // Carbon balance
   double preTotalC;
   double postTotalC;
-  double finalC;
   double inputsC;
   double outputsC;
 
   // Nitrogen balance
   double preTotalN;
   double postTotalN;
-  double finalN;
   double inputsN;
   double outputsN;
 
-  // Stock clamping adjustments: mass added by clamping negative stocks to zero.
-  // Positive values indicate that stocks were clamped (mass was "created" to
-  // prevent negative pools).
-  double clampedC;
-  double clampedN;
-
   // Checks
   double deltaC;
   double deltaN;
@@ -46,11 +43,10 @@ void updateBalanceTrackerPreUpdate(void);
 
 void updateBalanceTrackerPostUpdate(void);
 
-void updateBalanceTrackerPostClamp(void);
-
 void initBalanceTracker(void);
 
-// Returns 1 if a balance residual is nonfinite or has magnitude >= EPS
+// Returns 1 if a balance residual is nonfinite or has magnitude >=
+// BALANCE_TOLERANCE
 int balanceResidualFailed(double residual);
 
 // Exits with EXIT_CODE_INTERNAL_ERROR if the carbon or nitrogen balance fails

src/sipnet/sipnet.c

@@ -1796,15 +1796,12 @@ void updatePoolsAndBalance() {
   // correct for harvest termination
   checkForMortality();
 
-  // Calc total C and N after pool updates
-  updateBalanceTrackerPostUpdate();
-
-  // Verify none of our stocks have gone negative (set any that are to zero).
+  // Verify none of our stocks have gone negative; exits on a real violation,
+  // and zeroes values within rounding error of zero
   ensureNonNegativeStocks();
 
-  // Calc total C and N after non-negative check, and total system inputs and
-  // outputs
-  updateBalanceTrackerPostClamp();
+  // Calc total C and N after pool updates, and total system inputs and outputs
+  updateBalanceTrackerPostUpdate();
 
   // Perform balance check
   checkBalance();

tests/sipnet/test_modeling/testBalanceCheck.c

@@ -18,13 +18,11 @@ static void setResidual(Element element, double residual) {
   if (element == CARBON) {
     balanceTracker.preTotalC = 0.0;
     balanceTracker.postTotalC = residual;
-    balanceTracker.finalC = residual;
     balanceTracker.inputsC = 2.0;
     balanceTracker.outputsC = 2.0;
   } else {
     balanceTracker.preTotalN = 0.0;
     balanceTracker.postTotalN = residual;
-    balanceTracker.finalN = residual;
     balanceTracker.inputsN = 2.0;
     balanceTracker.outputsN = 2.0;
   }
@@ -83,18 +81,20 @@ static int expectPass(Element element, double residual, const char *label) {
 
 static int testElement(Element element) {
   int status = 0;
-  double justBelow = nextafter(EPS, 0.0);
+  double justBelow = nextafter(BALANCE_TOLERANCE, 0.0);
 
   // Balanced and below-threshold residuals are numerical noise
   status |= expectPass(element, 0.0, "zero");
-  status |= expectPass(element, 0.5 * EPS, "half-EPS positive");
-  status |= expectPass(element, -0.5 * EPS, "half-EPS negative");
-  status |= expectPass(element, justBelow, "just-below-EPS positive");
-  status |= expectPass(element, -justBelow, "just-below-EPS negative");
+  status |=
+      expectPass(element, 0.5 * BALANCE_TOLERANCE, "half-tolerance positive");
+  status |=
+      expectPass(element, -0.5 * BALANCE_TOLERANCE, "half-tolerance negative");
+  status |= expectPass(element, justBelow, "just-below-tolerance positive");
+  status |= expectPass(element, -justBelow, "just-below-tolerance negative");
 
   // At and above the threshold, in both directions
-  status |= expectFatal(element, EPS, "at-EPS positive");
-  status |= expectFatal(element, -EPS, "at-EPS negative");
+  status |= expectFatal(element, BALANCE_TOLERANCE, "at-tolerance positive");
+  status |= expectFatal(element, -BALANCE_TOLERANCE, "at-tolerance negative");
   status |= expectFatal(element, 1e-3, "large positive");
   status |= expectFatal(element, -1e-3, "large negative");
 
@@ -108,14 +108,14 @@ static int testElement(Element element) {
 
 static int testResidualPredicate(void) {
   int status = 0;
-  double justBelow = nextafter(EPS, 0.0);
+  double justBelow = nextafter(BALANCE_TOLERANCE, 0.0);
 
   status |= balanceResidualFailed(0.0) != 0;
   status |= balanceResidualFailed(-0.0) != 0;
   status |= balanceResidualFailed(justBelow) != 0;
   status |= balanceResidualFailed(-justBelow) != 0;
-  status |= balanceResidualFailed(EPS) != 1;
-  status |= balanceResidualFailed(-EPS) != 1;
+  status |= balanceResidualFailed(BALANCE_TOLERANCE) != 1;
+  status |= balanceResidualFailed(-BALANCE_TOLERANCE) != 1;
   status |= balanceResidualFailed(NAN) != 1;
   status |= balanceResidualFailed(INFINITY) != 1;
   status |= balanceResidualFailed(-INFINITY) != 1;
@@ -132,14 +132,11 @@ static int testBalancedStep(void) {
   initBalanceTracker();
   balanceTracker.preTotalC = 1000.0;
   balanceTracker.postTotalC = 1006.0;
-  balanceTracker.finalC = 1006.0;
   balanceTracker.inputsC = 10.0;
   balanceTracker.outputsC = 4.0;
   balanceTracker.preTotalN = 50.0;
-  balanceTracker.postTotalN = 49.5;
-  balanceTracker.finalN = 50.25;
-  balanceTracker.clampedN = 0.75;
-  balanceTracker.inputsN = 1.75;  // includes clamped
+  balanceTracker.postTotalN = 50.25;
+  balanceTracker.inputsN = 1.75;
   balanceTracker.outputsN = 1.5;
 
   really_exit = 0;

tests/sipnet/test_system/testCompleteHarvest.c

@@ -13,8 +13,6 @@ static ModelParams *testParams;
 static void balanced(void) {
   near(balanceTracker.deltaC, 0, "C balance");
   near(balanceTracker.deltaN, 0, "N balance");
-  near(balanceTracker.clampedC, 0, "C clamping");
-  near(balanceTracker.clampedN, 0, "N clamping");
 }
 
 static void processEvents(void) {

tests/sipnet/test_system/testLimitations.c

@@ -124,8 +124,6 @@ static void limitationCase(int harvest, int leafOn, int leafOff, int fertilizer,
   updateState();
   near(balanceTracker.deltaC, 0, "C balance");
   near(balanceTracker.deltaN, 0, "N balance");
-  near(balanceTracker.clampedC, 0, "C clamping");
-  near(balanceTracker.clampedN, 0, "N clamping");
   near(envi.plantLeafC, 10 - 2 * harvest + newLeaves - shedLeaves + growth,
        "plantLeafC");
   near(envi.litterC, 100 + 11 * harvest + shedLeaves + 3 * fertilizer,