Refactoring flux files

[seantrim]

Make sure all the relevant boxes are checked (and only check the box if you actually completed the step):

• [ ] Closes #xxx (identify the issue associated with this PR)
• [x] Code passes standard test cases (results are either bit-for-bit identical, or differences are explained in the PR comment)
• [ ] New tests added (describe which tests were performed to test the changes)
• [ ] Science test figures (add figures to PR comment and describe the tests)
• [ ] Checked that the new code conforms to the SUMMA coding conventions
• [x] Describe the change in the release notes (use either ./summa/docs/whats-new.md or ./summa/docs/minor-changes.md depending on what changed)

Refactoring of the flux related files vegNrgFlux.f90, ssdNrgFlux.f90, vegLiqFlux.f90, snowLiqFlux.f90, soilLiqFlux.f90, groundwatr.f90, bigAquifer.f90, and computFlux.f90 was performed. This includes cleanup, shortening, and updates to conform to the modern Fortran standard. Overall, the flux files have been shortened by nearly 900 lines. Updates to Fortran arithmetic expressions were also made to improve accuracy while reducing computational cost (such as using integer exponents and the intrinsic sqrt function where possible). There was also a bug with the initialization of C_r in vegNrgFlux.f90 that has been fixed (C_r was unintentionally assigned a literal of default real type -- it is now assigned a value of type rkind). Additionally, a script has been added to simplify the compilation process of SUMMA (summa/build/compile_SUMMA.sh), which was helpful when recompiling during the testing process.

The laugh tests were used to confirm functionality of the proposed changes. Apart from the Wigmosta tests, all laugh tests give bit identical results compared to the original develop branch. For the Wigmosta tests, the change in variable output does not exceed $10^{-7}$. This change is due to reduced round off error in several arithmetic expressions.

[wknoben]

Hi Sean,

Thanks for this. Three things:

1. Can you confirm that you ran the full suite of SUMMA test cases (https://www.hydroshare.org/resource/e4df3065264341a48017f6d2f4ad030c/) and not just the laugh tests (https://github.com/CH-Earth/laughTests)? The full suite covers more of the model behavior than the laugh tests do.

2. Can you add an example (specific lines that were changed and a plot of a relevant flux or state variable) for the Wigmosta case? I don't doubt that you're right but it would be good to have an accurate record of which changes in the code caused these changes in output in case we ever need to get back to it. Also, how do we know round off error is reduced in the new implementation?

3. W.r.t the compile script, does that need any additional comments to indicate (1) to novice users that this is not fully automatic but may need manual updates to library paths, and (2) what the prerequisites for use are?

[seantrim]

Hi Wouter,

Thanks for your helpful comments. To respond to your points:

1) Originally I only ran the laugh tests. However, I have expanded things to the full test suite which includes the laugh tests (synthetic tests) and the tests based on field data (i.e., all 35 tests from https://www.hydroshare.org/resource/e4df3065264341a48017f6d2f4ad030c/). A summary of differences between output from the develop and develop_refactor branches can be found here. The file lists the maximum absolute and percentage differences detected for each variable for all 35 tests (this was obtained using a modified version of the check_bit_4_bit.py script). For percentage differences, it is possible to get division by zero if the variable changed from a value of zero (corresponding to values of inf or nan). However, the absolute differences are small for those cases. Plots of variables that showed relatively larger differences are available here (this includes the Wigmosta tests which are tests 11 and 12). Overall model behaviour matched very closely between the develop and develop refactor branches.

2) The differences in output were caused by updates to Fortran arithmetic in vegNrgFlux.f90 in multiple spots. For example, in line 3178 in the develop branch we have: dStabilityCorrection_dRich = (-16._rkind) * 0.5_rkind*(1._rkind - 16._rkind*RiBulk)**(-0.5_rkind) which has been updated to dStabilityCorrection_dRich = -8._rkind/sqrt(1._rkind - 16._rkind*RiBulk) In the above example, there is a multiplication (-16._rkind) * 0.5_rkind = -8._rkind, and a square root (1._rkind - 16._rkind*RiBulk)**(-0.5_rkind) = 1/sqrt(1._rkind - 16._rkind*RiBulk) that have been simplified. Eliminating the multiplication reduces the number of floating point operations, reducing both computational expense and the accumulation of round off error. Using the sqrt Fortran intrinsic is faster than raising to the power of 1/2 because the intrinsic applies a special algorithm that exists for square roots (see https://stackoverflow.com/questions/29098168/why-does-0-5-appear-to-be-more-efficient-than-sqrt for a speed test). Also, the sqrt intrinsic avoids the storage limitations on the 0.5_rkind literal which may have a small error in the least significant digit due to finite precision. Accordingly, the result using the sqrt intrinsic will in general be slightly more accurate (and faster) than the exponent operation. Related to this issue is the implementation of integer valued exponents. For instance, x**3 is more accurate than x**3.0_rkind for computing $x^3$ because the integer literal is stored exactly whereas 3.0_rkind is limited to finite precision (it will be 3 to within round off error). Also, x**3 will be interpreted as three multiplications which is less expensive than x**3.0_rkind (interpreted as exp(3.0_rkind*log(x))). A related discussion can be found at https://www.personal.psu.edu/jhm/f90/lectures/7.html. Integer exponents were implemented in the develop refactor branch where possible.

3) I've added additional comments on prerequisites and usage in the compile_SUMMA.sh script.

[wknoben]

Thanks for the updates.

1. I noticed that the discussion on sqrt() vs **0.5 refers to R instead of Fortran, though playing around with a test function in an online Fortran environment seems to confirm that sqrt() is considerably faster than **0.5 (3 tests each, execution time 7.70, 7.86, 7.79 for 0.5; 2.14, 2.13, 2.21 for sqrt()).

Code

``` Program test_sqrt_methods implicit none integer :: N,i real :: T1,T2,x,y ! initiate time call cpu_time(T1) ! Run square root N = 1000000000 x = 10 do i = 1,N !y = x**(0.5) y = sqrt(x) end do ! End time call cpu_time(T2) print *, T1,T2, T2-T1 End Program test_sqrt_methods ``` Tested here: https://www.onlinegdb.com/online_fortran_compiler

2. I'll have a look at the code changes and leave specific comments if I see anything. Given what kind of work this is it'd be good if @martynpclark can approve these too.

3. Last, for future work it might make sense to separate cosmetic changes and code changes into separate PRs. That'll make it easier to see those changes that fundamentally affect model simulations.

[andywood]

I think it's fine to retain some old comments and testing blocks unless something is absolutely and cleanly finalized and never needs to be seen again. We're not developing software for IBM, but rather hydrology models, and it's not terrible for the codes to keep some of the personality and thought process of their creators & developers. E.g., at some point someone wiped out Martyn's 'dirs_und_philes' subroutine in place of 'dirs_and_files' ... was that really necessary? Good to clean up code, fix bugs, improve documentation, remove clutter ... but I'm not sure there's a need for complete sanitization.

The rkind/i4b issue -- there was a search/replace when switching the model from fixed types to rkind, and it's reasonable that there are places where a fixed type makes sense. indexes are one of those cases (versus state variables). probably needs a case by case review when inconsistencies are found.

On Tue, May 30, 2023 at 4:32 PM Wouter Knoben @.***> wrote:

@.**** commented on this pull request.

In build/source/engine/soilLiqFlx.f90 https://github.com/CH-Earth/summa/pull/525#discussion_r1210888516:

• ! (test derivative)
• if(testDeriv)then
• xConst = LH_fus/(gravity*Tfreeze) ! m K-1 (NOTE: J = kg m2 s-2)
• vTheta = scalarVolFracIceTrial + scalarVolFracLiqTrial
• volLiq = volFracLiq(xConst*(scalarTempTrial+dx - Tfreeze),vGn_alpha,theta_res,theta_sat,vGn_n,vGn_m)
• volIce = vTheta - volLiq
• effSat = (volLiq - theta_res)/(theta_sat - volIce - theta_res)
• psiLiq = matricHead(effSat,vGn_alpha,0._rkind,1._rkind,vGn_n,vGn_m) ! use effective saturation, so theta_res=0 and theta_sat=1
• hydCon = hydCond_psi(psiLiq,scalarSatHydCond,vGn_alpha,vGn_n,vGn_m)
• call iceImpede(volIce,f_impede,iceImpedeFac,dIceImpede_dLiq)
• hydIce = hydCon*iceImpedeFac
• print*, 'test derivative: ', (psiLiq - scalarMatricHeadTrial)/dx, dPsiLiq_dTemp
• print*, 'test derivative: ', (hydCon - hydCond_noIce)/dx, dHydCondMicro_dTemp
• print*, 'test derivative: ', (hydIce - scalarHydCond)/dx, dHydCond_dTemp
• print, 'press any key to continue'; read(,*) ! (alternative to the deprecated 'pause' statement)
• end if ! testing the derivative

More old testing code. Worth keeping? @martynpclark https://github.com/martynpclark

— Reply to this email directly, view it on GitHub https://github.com/CH-Earth/summa/pull/525#pullrequestreview-1451922950, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABIKARIACJSJZWRTTON3HYLXIZYQLANCNFSM6AAAAAAYNDYFHE . You are receiving this because you are subscribed to this thread.Message ID: @.***>

[wknoben]

Perfect. For posterity, can you add why this can be removed (e.g., code moved to function X, functionality deprecated)?

[seantrim]

Thanks Wouter for looking things over! That exfiltration check in groundwatr.f90 was previously used for debugging purposes but is no longer needed (confirmed with @martynpclark during a meeting). I have amended the whats-new.md file with a statement explaining that several print statements previously used for debugging were removed from the flux files because they were no longer required.

[wknoben]

Thanks for the updates @seantrim!

@martynpclark @andywood I'll leave this open for a few more days to give you a chance to comment if you want. If not I'll merge this some time next week.