From c3d5b477ac85e9b27e099b75a1a733eadd165489 Mon Sep 17 00:00:00 2001 From: priyanshusky0 Date: Wed, 28 Oct 2026 10:21:52 +0530 Subject: [PATCH 1/2] Remove clamping from balance checks (#416) ensureNonNegative now exits on real violations, so the balance tracker no longer needs to account for mass added by clamping. Drop the clamped and final totals, take the post-update totals after ensureNonNegativeStocks(), and compute system inputs and outputs in updateBalanceTrackerPostUpdate(). Widen the balance tolerance to 10 * EPS to allow for several sub-EPS stocks being zeroed in the same step. --- src/sipnet/balance.c | 75 +++++-------------- src/sipnet/balance.h | 22 +++--- src/sipnet/sipnet.c | 11 +-- tests/sipnet/test_modeling/testBalanceCheck.c | 31 ++++---- .../sipnet/test_system/testCompleteHarvest.c | 2 - tests/sipnet/test_system/testLimitations.c | 2 - 6 files changed, 44 insertions(+), 99 deletions(-) diff --git a/src/sipnet/balance.c b/src/sipnet/balance.c index f71775c6..8c326038 100644 --- a/src/sipnet/balance.c +++ b/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,12 +64,6 @@ 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) { @@ -103,12 +71,8 @@ void initBalanceTracker(void) { // 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,48 +86,45 @@ 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; } @@ -171,16 +132,14 @@ void checkBalance(void) { 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"); diff --git a/src/sipnet/balance.h b/src/sipnet/balance.h index 3a8723d9..9afb97dd 100644 --- a/src/sipnet/balance.h +++ b/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 diff --git a/src/sipnet/sipnet.c b/src/sipnet/sipnet.c index 1ba894fa..6a1a6e8b 100644 --- a/src/sipnet/sipnet.c +++ b/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(); diff --git a/tests/sipnet/test_modeling/testBalanceCheck.c b/tests/sipnet/test_modeling/testBalanceCheck.c index d46b80ef..c2ec6d3a 100644 --- a/tests/sipnet/test_modeling/testBalanceCheck.c +++ b/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; diff --git a/tests/sipnet/test_system/testCompleteHarvest.c b/tests/sipnet/test_system/testCompleteHarvest.c index d60c6d5a..e90544b6 100644 --- a/tests/sipnet/test_system/testCompleteHarvest.c +++ b/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) { diff --git a/tests/sipnet/test_system/testLimitations.c b/tests/sipnet/test_system/testLimitations.c index 4525c2ec..0f0a3b93 100644 --- a/tests/sipnet/test_system/testLimitations.c +++ b/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, From 5cf89ab88a481a1252ec82f56a87aa51b3f579f2 Mon Sep 17 00:00:00 2001 From: priyanshusky0 Date: Wed, 28 Oct 2026 10:59:05 +0530 Subject: [PATCH 2/2] Add CHANGELOG entry for #425 --- docs/CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 57a4a002..5ef8fda3 100644 --- a/docs/CHANGELOG.md +++ b/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