Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

fix(UzfCellGroup): SQUARE_GWET option records ET as inflow instead of outflow #1951

Open
wants to merge 24 commits into
base: develop
Choose a base branch
from

Conversation

emorway-usgs
Copy link
Contributor

GWET values calculated using the SQUARE_GWET option in the UZF options block had the wrong sign.

  • Replaced section above with description of pull request
  • Resolves discussion UZF-GWET term on volume budget #1904
  • Adds a new test
  • Ran ruff on new and modified python scripts in autotests subdirectory.
  • Formatted new and modified Fortran source files with fprettify
  • Updated develop.tex with a plain-language description of the bug fix, change, feature; required for changes that may affect users

@emorway-usgs
Copy link
Contributor Author

Some further notes on GWET as calculated by UZF:

For the SQUARE_GWET option in UZF, thcof (shorthand for "total uzf hcof contribution to GWF model") in the etfunc_nlin is not a function of head except over the smoothing interval at the extinction depth. This nonlinearity is handled directly through the Jacobian in the uzf8 source file:

    ! -- store derivative for Newton addition to equations in _fn()
    this%deriv(i) = uzderiv + derivgwet

Also, note that in the call to etfunc_lin, thcof is updated and returned to simgwet. However, in the function etfunc_nlin, thcof is not updated and remains zero. As result, in the following code, found in simgwet

if (igwetflag == 1) then
  et = etfunc_lin(s, x, c, det, trhs, thcof, hgwf, &
                  this%celtop(icell), this%celbot(icell))
else if (igwetflag == 2) then
  et = etfunc_nlin(s, x, c, det, trhs, thcof, hgwf)
  signflip = -1
end if
! this%gwet(icell) = et * this%uzfarea(icell)
trhs = trhs * this%uzfarea(icell)
thcof = thcof * this%uzfarea(icell)
this%gwet(icell) = (trhs - (thcof * hgwf))

the last line will always be trhs (shorthand for "total uzf rhs contribution to GWF model") minus 0.0 since thcof will remain at its initialized value of DZERO for the SQUARE_GWET option. As a result, %gwet will be positive when calculated using the SQUARE_GWET option but will be negative when the LINEAR_GWET option is used.

Copy link
Contributor

@langevin-usgs langevin-usgs left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nice job digging this out. Couple of comments included in this review.

doc/ReleaseNotes/develop.tex Outdated Show resolved Hide resolved
@@ -875,18 +877,19 @@ subroutine simgwet(this, igwetflag, icell, hgwf, trhs, thcof, det)
this%celtop(icell), this%celbot(icell))
else if (igwetflag == 2) then
et = etfunc_nlin(s, x, c, det, trhs, thcof, hgwf)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there any way to adjust etfunc_lin and etfunc_nlin so that they return hcof, rhs, and rate in a way that is consistent with one another? It's unfortunate that there has to be sign shenanigans on the outside of the calls to these routines. Just a thought.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Have shifted etfunc_lin and etfunc_nlin over to a dedicated file that is included in the unit testing.

@emorway-usgs emorway-usgs force-pushed the fix_square_gwet branch 3 times, most recently from 27b869e to b9e7a00 Compare July 23, 2024 15:45
@emorway-usgs emorway-usgs force-pushed the fix_square_gwet branch 3 times, most recently from 08f6c61 to 32257db Compare August 1, 2024 00:01
@emorway-usgs emorway-usgs force-pushed the fix_square_gwet branch 2 times, most recently from 5026191 to 5a61c2b Compare August 8, 2024 17:05
@emorway-usgs emorway-usgs force-pushed the fix_square_gwet branch 2 times, most recently from a7ba206 to cba0802 Compare August 15, 2024 21:44
@emorway-usgs emorway-usgs marked this pull request as draft September 4, 2024 14:20
@emorway-usgs emorway-usgs force-pushed the fix_square_gwet branch 3 times, most recently from 4d6e43b to eca2a78 Compare October 23, 2024 02:21
@emorway-usgs emorway-usgs force-pushed the fix_square_gwet branch 2 times, most recently from e03ff70 to 6a306d5 Compare November 5, 2024 18:02
@emorway-usgs emorway-usgs marked this pull request as ready for review November 6, 2024 15:48
@emorway-usgs
Copy link
Contributor Author

Added a new test problem with checks to ensure the GWET is being removed and that the resulting heads are as expected. The test compares results when using the classic LINEAR_GWET vs SQUARE_GWET

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants