Ptl mbal testing #21

peterlafollette · 2024-03-25T21:24:23Z

[Short description explaining the high-level reason for the pull request]

Additions

After the calibration workstream experienced some errors with LASAM, we decided it was a good idea to run LASAM with a large number of parameter sets and debug all errors that arose. In this branch, I've tested LASAM with 42k randomly generated parameter sets, including humid and arid forcing data, in the framework and in standalone mode, and in one and 3 layer scenarios. Currently the model runs successfully with each parameter set.

Debugging mostly had to do with implementing many small changes to address unlikely scenarios that we did not anticipate earlier in the model's development. While the bulk of the code changes address such errors (brief comments are given describing why the code was added in most cases), this push also has the implementation of an adaptive integral to calculate the capillary drive G, as well as a redefinition of field capacity as a head rather than as a fixed relative water content. Both of these changes were developed as part of LGARTO model development.

Removals

Removals focus on lines of code that yielded errors in some tested cases.

Changes

Testing

Testing was achieved by running the model with 42k unique parameter sets.

Screenshots

Notes

Todos

Checklist

Testing checklist

Target Environment support

Windows
Linux
Browser

Accessibility

Keyboard friendly
Screen reader friendly

Other

peterlafollette · 2024-03-25T21:30:18Z

I just realized, the unit tests need to be updated because the calculation of capillary drive G is now slightly different ...

peterlafollette · 2024-03-25T22:06:49Z

remaking PR with updated unit test.

peterlafollette · 2024-03-25T22:13:12Z

Remade PR with unit test updated. I have not yet updated the synthetic tests and will likely suggest a new (simpler) way forward for unit testing in a meeting tomorrow. Essentially, I think we should just do unit testing on the single, simple test, and then provide a jupyter notebook that allows for more widespread testing.

…ance errors and infinite loop errors that became apparent when LASAM was run with many parameter sets. Also updated field capacity formulation to use a head rather than a relative water content, and updated caluclation for capillary drive G to use an adaptive integral, both as in LGARTO. Rebasing to incorporate recent changes.

…ance errors and infinite loop errors that became apparent when LASAM was run with many parameter sets. Also updated field capacity formulation to use a head rather than a relative water content, and updated caluclation for capillary drive G to use an adaptive integral, both as in LGARTO. Cleaned up some commenting and redundancy after edits were made.

src/lgar.cxx

ajkhattak · 2024-04-02T16:44:13Z

src/lgar.cxx

+  int iter = 0;
+	while (fabs(mass_balance_error - tolerance) > 1.E-10) {
+    iter++;
+    if (iter>10000000) {


this is a big number, how often does it happen that the code does not converge? Instead of aborting we need to fix this (but not in this PR)....

In its current state, this code never triggers. I used it while debugging, but kept it in because it seems good to not let the loop get to this many iterations.

I believe the code should converge in 100s of iterations (<1000 I would guess). Is it (the number 10^6) related to avoiding an infinite loop? If so, we should use a smaller number (say, 1000-5000). If the convergence does not happen in 1000-5000 iterations, we should fix it in a better way. With 10^6, the code might converge in 100,000 iterations -and we will not know about it -- which I don't think is acceptable.

Good point. Rewrote such that the loop terminates after 1e4 iterations (rather than the original 1e7) and puts the remaining mass balance error into AET, rather than aborting, similar to how we do it in the other mass balance loop (in the function lgar_fix_dry_over_wet_wetting_fronts).

src/lgar.cxx

src/soil_funcs.cxx

ajkhattak · 2024-04-04T15:40:01Z

src/lgar.cxx

@@ -1587,11 +1643,15 @@ extern void lgar_wetting_fronts_cross_layer_boundary(int num_layers,
      current->dzdt_cm_per_h = 0;
      current->to_bottom = TRUE;
      next->to_bottom = FALSE;
+
+      if (isinf(next->depth_cm)){//in some rare cases, due to (very) different shapes of soil water rentention curves in adjacent soil layers near saturation, the WF that crossed a layer bdy can yield a theta value identical to the one below it, resulting in division by 0 and an infinite WF depth. 


the code should never assign infinity to a wetting front depth, instead of checking it here, the code should handle it early on. I think the right place is:

double mbal_Z_correction = mbal_correction / (theta_new - next_to_next->theta);

We should also add an assertion
assert (wetting_front->depth_cm > domain_depth)

I updated the code as you suggested, fixing the error at mbal_Z_correction rather than when we assign the WF depth, but I don't think we should add that assertion. This is because we have a whole function for crossing the domain lower boundary that handles that case.

okay, so I guess this was the only place where WF depth was crossing the domain boundary depth, correct?

Well, generally WF depth can cross the domain boundary depth any time WF depth is updated, so this is not the only place. However, the cases of wetting fronts crossing a layer boundary, wetting fronts crossing a domain boundary, and wetting fronts merging with one another are handled specially by an individual function for each of those, and in specific places. The code you commented on here is in the function where WFs cross a layer boundary, so in this function we don't deal with the case where a WF crosses the domain bdy.

I should add, this is probably the only place where WF depth can be calculated as infinite. This is because this is probably the only place where we have a delta theta in the denominator used to update WF depth (aside from in dzdt function, where we also deal with the issue of 0 in the denominator).

ajkhattak · 2024-04-04T15:49:14Z

src/lgar.cxx

+
+    //so the above code that avoids infinite loops by breaking in certain cases works, but a fringe case must be addressed. 
+    //When theta is close to theta_r, the loop "while (delta_mass > tolerance)" will attempt to increase psi until mass balance closure occurs, but if the 
+    //AET value is large enough to make theta less than theta_r (or if dzdt was large enough to make the WF move such that the updated theta would be less than theta_r),
+    //the loop would never be able to achieve mass balance closure, and it will break when the one of above conditions is met.


please remove the comments that are no longer relevant such as theta less than theta_r, we don't allow this condition anymore, correct?

You're correct; we don't allow this condition. This comment is relevant however, because it explains that the following code addresses the case where theta would go below theta_r and instead prevents this.

ajkhattak · 2024-04-04T15:56:31Z

src/lgar.cxx

@@ -2432,6 +2537,17 @@ extern double lgar_theta_mass_balance(int layer_num, int soil_num, double psi_cm
    if (psi_cm_loc <= 0 && psi_cm_loc_prev < 1E-50) break;


1.0E-50 is 0.0

I agree that 1e-50 is basically 0. I'm testing with 0 instead, I think that it should be fine.

… less than theta_e

…ll never yield theta<theta_r

src/lgar.cxx

ajkhattak

I am sure with these changes the model is more stable now. Thanks for working on it.

peterlafollette requested a review from ajkhattak March 25, 2024 21:24

peterlafollette closed this Mar 25, 2024

peterlafollette reopened this Mar 25, 2024

Peter La Follette added 3 commits April 2, 2024 10:04

updated unit test

48b0f56

ajkhattak requested changes Apr 2, 2024

View reviewed changes

peterlafollette force-pushed the PTL_mbal_testing branch from 933a66f to 48b0f56 Compare April 2, 2024 18:23

updating with suggestions from PR review

032072d

ajkhattak reviewed Apr 4, 2024

View reviewed changes

Peter La Follette added 4 commits April 4, 2024 10:04

Updated lgar.cxx in response to a few comments

afc7673

Removed a few comments that were not necessary

ed41fc5

when theta updates, it is now enforced to be greater than theta_r and…

5cc3d70

… less than theta_e

increased clarity of comment detailing how lgar_theta_mass_balance wi…

c2f0c79

…ll never yield theta<theta_r

ajkhattak reviewed Apr 4, 2024

View reviewed changes

src/lgar.cxx Show resolved Hide resolved

ajkhattak approved these changes Apr 4, 2024

View reviewed changes

ajkhattak merged commit 4c21c46 into master Apr 4, 2024
4 checks passed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Ptl mbal testing #21

Ptl mbal testing #21

peterlafollette commented Mar 25, 2024 •

edited

Loading

peterlafollette commented Mar 25, 2024

peterlafollette commented Mar 25, 2024

peterlafollette commented Mar 25, 2024

ajkhattak Apr 2, 2024

peterlafollette Apr 2, 2024

ajkhattak Apr 4, 2024

peterlafollette Apr 4, 2024 •

edited

Loading

ajkhattak Apr 4, 2024

peterlafollette Apr 4, 2024 •

edited

Loading

ajkhattak Apr 4, 2024

peterlafollette Apr 4, 2024

peterlafollette Apr 4, 2024

ajkhattak Apr 4, 2024

peterlafollette Apr 4, 2024

ajkhattak Apr 4, 2024

peterlafollette Apr 4, 2024 •

edited

Loading

ajkhattak left a comment

		@@ -2432,6 +2537,17 @@ extern double lgar_theta_mass_balance(int layer_num, int soil_num, double psi_cm
		if (psi_cm_loc <= 0 && psi_cm_loc_prev < 1E-50) break;

Ptl mbal testing #21

Ptl mbal testing #21

Conversation

peterlafollette commented Mar 25, 2024 • edited Loading

Additions

Removals

Removals focus on lines of code that yielded errors in some tested cases.

Changes

Testing

Screenshots

Notes

Todos

Checklist

Testing checklist

Target Environment support

Accessibility

Other

peterlafollette commented Mar 25, 2024

peterlafollette commented Mar 25, 2024

peterlafollette commented Mar 25, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

peterlafollette Apr 4, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

peterlafollette Apr 4, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

peterlafollette Apr 4, 2024 • edited Loading

Choose a reason for hiding this comment

ajkhattak left a comment

Choose a reason for hiding this comment

peterlafollette commented Mar 25, 2024 •

edited

Loading

peterlafollette Apr 4, 2024 •

edited

Loading

peterlafollette Apr 4, 2024 •

edited

Loading

peterlafollette Apr 4, 2024 •

edited

Loading