diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index 744e43ec16..df0450016b 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -38,6 +38,7 @@ module FATESPlantRespPhotosynthMod use FatesInterfaceTypesMod, only : hlm_parteh_mode use FatesInterfaceTypesMod, only : numpft use FatesInterfaceTypesMod, only : nleafage + use FatesUtilsMod, only : QuadraticRoots => QuadraticRootsSridharachary use EDParamsMod, only : maxpft use EDParamsMod, only : nlevleaf use EDParamsMod, only : nclmax @@ -1379,7 +1380,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in aquad = theta_psii bquad = -(qabs + jmax) cquad = qabs * jmax - call quadratic_f (aquad, bquad, cquad, r1, r2) + call QuadraticRoots(aquad, bquad, cquad, r1, r2) je = min(r1,r2) ! Initialize intercellular co2 @@ -1409,7 +1410,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in aquad = theta_cj_c3 bquad = -(ac + aj) cquad = ac * aj - call quadratic_f (aquad, bquad, cquad, r1, r2) + call QuadraticRoots(aquad, bquad, cquad, r1, r2) agross = min(r1,r2) else @@ -1440,13 +1441,13 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in aquad = theta_cj_c4 bquad = -(ac + aj) cquad = ac * aj - call quadratic_f (aquad, bquad, cquad, r1, r2) + call QuadraticRoots(aquad, bquad, cquad, r1, r2) ai = min(r1,r2) aquad = theta_ip bquad = -(ai + ap) cquad = ai * ap - call quadratic_f (aquad, bquad, cquad, r1, r2) + call QuadraticRoots(aquad, bquad, cquad, r1, r2) agross = min(r1,r2) end if @@ -1485,7 +1486,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in (2.0*stomatal_intercept_btran + term * & (1.0 - medlyn_slope(ft)* medlyn_slope(ft) / vpd)) * term - call quadratic_f (aquad, bquad, cquad, r1, r2) + call QuadraticRoots(aquad, bquad, cquad, r1, r2) gs_mol = max(r1,r2) else if ( stomatal_model == ballberry_model ) then !stomatal conductance calculated from Ball et al. (1987) @@ -1494,7 +1495,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in cquad = -gb_mol*(leaf_co2_ppress*stomatal_intercept_btran + & bb_slope(ft)*anet*can_press * ceair/ veg_esat ) - call quadratic_f (aquad, bquad, cquad, r1, r2) + call QuadraticRoots(aquad, bquad, cquad, r1, r2) gs_mol = max(r1,r2) end if @@ -1922,145 +1923,6 @@ end function fth25_f ! ===================================================================================== - subroutine quadratic_f (a, b, c, r1, r2) - ! - ! !DESCRIPTION: - !==============================================================================! - !----------------- Solve quadratic equation for its two roots -----------------! - !==============================================================================! - ! This solution is mostly derived from: - ! Press WH, Teukolsky SA, Vetterling WT, Flannery BP. 1992. Numerical Recipes - ! in Fortran77: The Art of Scientific Computing. 2nd edn. Cambridge - ! University Press, Cambridge UK, ISBN 0-521-43064-X. - ! Available at: http://numerical.recipes/oldverswitcher.html, section 5.6. - ! - ! !REVISION HISTORY: - ! 4/5/10: Adapted from /home/bonan/ecm/psn/An_gs_iterative.f90 by Keith Oleson - ! 7/23/16: Copied over from CLM by Ryan Knox - ! 12/30/23: Instead of issuing errors when a=0, solve the trivial cases too. - ! Check determinant sign, and stop the run when it is negative. - ! - ! !USES: - ! - ! !ARGUMENTS: - real(r8), intent(in) :: a,b,c ! Terms for quadratic equation - real(r8), intent(out) :: r1,r2 ! Roots of quadratic equation - ! - ! !LOCAL VARIABLES: - real(r8) :: discriminant ! Discriminant - real(r8) :: q ! Temporary term for quadratic solution - logical :: a_offzero ! Is a close to zero? - logical :: b_offzero ! Is b close to zero? - logical :: c_offzero ! Is c close to zero? - ! ! Local constants: - real(r8), parameter :: discard = 1.e36_r8 ! Large number for discarding answer - !------------------------------------------------------------------------------ - - ! Save logical tests. - a_offzero = abs(a) > nearzero - b_offzero = abs(b) > nearzero - c_offzero = abs(c) > nearzero - - if (a_offzero .and. ( b_offzero .or. c_offzero ) ) then - ! Quadratic equation with two non-zero solutions (but may be complex solutions) - discriminant = b*b - 4._r8 * a * c - - ! Proceed only when the discriminant is non-negative or only tiny negative - if (discriminant >= - nearzero) then - ! Coerce discriminant to non-negative - discriminant = max(0._r8,discriminant) - - ! Find q as in the numerical recipes. If b or c are non-zero, q cannot - ! be zero, no need for additional checks. - q = - 0.5_r8 * (b + sign(sqrt(discriminant),b)) - r1 = q / a - r2 = c / q - else - ! Negative discriminant, stop the run. - write (fates_log(),'(a)') '---~---' - write (fates_log(),'(a)') ' Fatal error!' - write (fates_log(),'(a)') ' Quadratic equation discriminant is negative.' - write (fates_log(),'(a)') '---~---' - write (fates_log(),'(a,1x,es12.5)') ' a = ',a - write (fates_log(),'(a,1x,es12.5)') ' b = ',b - write (fates_log(),'(a,1x,es12.5)') ' c = ',c - write (fates_log(),'(a,1x,es12.5)') ' discriminant = ',discriminant - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - else if (a_offzero) then - ! b and c are nearly zero. Both roots must be zero. - r1 = 0._r8 - r2 = 0._r8 - else if (b_offzero) then - ! "a" is not zero, not a true quadratic equation. Single root. - r1 = - c / b - r2 = discard - else - ! Both a and b are zero, this really doesn't make any sense and should never - ! happen. If it does, issue an error and stop the run. - write (fates_log(),'(a)') '---~---' - write (fates_log(),'(a)') ' Fatal error!' - write (fates_log(),'(a)') ' This solver expects ''a'' and/or ''b'' to be non-zero.' - write (fates_log(),'(a)') '---~---' - write (fates_log(),'(a,1x,es12.5)') ' a = ',a - write (fates_log(),'(a,1x,es12.5)') ' b = ',b - write (fates_log(),'(a,1x,es12.5)') ' c = ',c - write (fates_log(),'(a)') '---~---' - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - - return - end subroutine quadratic_f - - ! ==================================================================================== - - subroutine quadratic_fast (a, b, c, r1, r2) - ! - ! !DESCRIPTION: - !==============================================================================! - !----------------- Solve quadratic equation for its two roots -----------------! - ! THIS METHOD SIMPLY REMOVES THE DIV0 CHECK AND ERROR REPORTING ! - !==============================================================================! - ! Solution from Press et al (1986) Numerical Recipes: The Art of Scientific - ! Computing (Cambridge University Press, Cambridge), pp. 145. - ! - ! !REVISION HISTORY: - ! 4/5/10: Adapted from /home/bonan/ecm/psn/An_gs_iterative.f90 by Keith Oleson - ! 7/23/16: Copied over from CLM by Ryan Knox - ! - ! !USES: - ! - ! !ARGUMENTS: - real(r8), intent(in) :: a,b,c ! Terms for quadratic equation - real(r8), intent(out) :: r1,r2 ! Roots of quadratic equation - ! - ! !LOCAL VARIABLES: - real(r8) :: q ! Temporary term for quadratic solution - !------------------------------------------------------------------------------ - - ! if (a == 0._r8) then - ! write (fates_log(),*) 'Quadratic solution error: a = ',a - ! call endrun(msg=errMsg(sourcefile, __LINE__)) - ! end if - - if (b >= 0._r8) then - q = -0.5_r8 * (b + sqrt(b*b - 4._r8*a*c)) - else - q = -0.5_r8 * (b - sqrt(b*b - 4._r8*a*c)) - end if - - r1 = q / a - ! if (q /= 0._r8) then - r2 = c / q - ! else - ! r2 = 1.e36_r8 - ! end if - - end subroutine quadratic_fast - - - ! ==================================================================================== - subroutine UpdateCanopyNCanNRadPresent(currentPatch) ! --------------------------------------------------------------------------------- diff --git a/functional_unit_testing/math_utils/MathUtilsDriver.py b/functional_unit_testing/math_utils/MathUtilsDriver.py new file mode 100644 index 0000000000..4add288126 --- /dev/null +++ b/functional_unit_testing/math_utils/MathUtilsDriver.py @@ -0,0 +1,172 @@ +# ======================================================================================= +# +# For usage: $python HydroUTestDriver.py --help +# +# This script runs unit tests on the hydraulics functions. +# +# +# ======================================================================================= + +import matplotlib as mpl +#mpl.use('Agg') +import matplotlib.pyplot as plt +from datetime import datetime +import argparse +#from matplotlib.backends.backend_pdf import PdfPages +import platform +import numpy as np +import os +import sys +import getopt +import code # For development: code.interact(local=dict(globals(), **locals())) +import time +import imp +import ctypes +from ctypes import * +from operator import add + + +#CDLParse = imp.load_source('CDLParse','../shared/py_src/CDLParse.py') +#F90ParamParse = imp.load_source('F90ParamParse','../shared/py_src/F90ParamParse.py') +PyF90Utils = imp.load_source('PyF90Utils','../shared/py_src/PyF90Utils.py') + + +#from CDLParse import CDLParseDims, CDLParseParam, cdl_param_type +#from F90ParamParse import f90_param_type, GetSymbolUsage, GetPFTParmFileSymbols, MakeListUnique + +from PyF90Utils import c8, ci, cchar, c8_arr, ci_arr + +# Load the fortran objects via CTYPES + +f90_unitwrap_obj = ctypes.CDLL('bld/UnitWrapMod.o',mode=ctypes.RTLD_GLOBAL) +f90_constants_obj = ctypes.CDLL('bld/FatesConstantsMod.o',mode=ctypes.RTLD_GLOBAL) +f90_fatesutils_obj = ctypes.CDLL('bld/FatesUtilsMod.o',mode=ctypes.RTLD_GLOBAL) + +# Alias the F90 functions, specify the return type +# ----------------------------------------------------------------------------------- + +neighbor_dist = f90_fatesutils_obj.__fatesutilsmod_MOD_getneighbordistance +#quadratic_f = f90_fatesutils_obj.__fatesutilsmod_MOD_quadratic_f +quadratic_roots = f90_fatesutils_obj.__fatesutilsmod_MOD_quadraticroots +quadratic_sroots = f90_fatesutils_obj.__fatesutilsmod_MOD_quadraticrootssridharachary + +# Some constants +rwcft = [1.0,0.958,0.958,0.958] +rwccap = [1.0,0.947,0.947,0.947] +pm_leaf = 1 +pm_stem = 2 +pm_troot = 3 +pm_aroot = 4 +pm_rhiz = 5 + +# These parameters are matched with the indices in FATES-HYDRO +vg_type = 1 +cch_type = 2 +tfs_type = 3 + +isoil1 = 0 # Top soil layer parameters (@BCI) +isoil2 = 1 # Bottom soil layer parameters + +# Constants for rhizosphere +watsat = [0.567, 0.444] +sucsat = [159.659, 256.094] +bsw = [6.408, 9.27] + +unconstrained = True + + +# ======================================================================================== +# ======================================================================================== +# Main +# ======================================================================================== +# ======================================================================================== + +def main(argv): + + # First check to make sure python 2.7 is being used + version = platform.python_version() + verlist = version.split('.') + + #if( not ((verlist[0] == '2') & (verlist[1] == '7') & (int(verlist[2])>=15) ) ): + # print("The PARTEH driver mus be run with python 2.7") + # print(" with tertiary version >=15.") + # print(" your version is {}".format(version)) + # print(" exiting...") + # sys.exit(2) + + # Read in the arguments + # ======================================================================================= + + # parser = argparse.ArgumentParser(description='Parse command line arguments to this script.') + # parser.add_argument('--cdl-file', dest='cdlfile', type=str, \ + # help="Input CDL filename. Required.", required=True) + # args = parser.parse_args() + + # Set number of analysis points + + # y = ax2 + bx + c + + a = [1,1,5,1.5] + b = [-2,7,10,3.2] + c = [1,12,3,1.1] + + cd_r1 = c_double(-9.0) + cd_r2 = c_double(-9.0) + + r1 = np.zeros([3,1]) + r2 = np.zeros([3,1]) + + for ic in range(len(a)): + + #iret = quadratic_f(c8(a[ic]),c8(b[ic]),c8(c[ic]),byref(cd_r1),byref(cd_r2)) + #r1[0] = cd_r1.value + #r2[0] = cd_r2.value + + iret = quadratic_roots(c8(a[ic]),c8(b[ic]),c8(c[ic]),byref(cd_r1),byref(cd_r2)) + r1[1] = cd_r1.value + r2[1] = cd_r2.value + + iret = quadratic_sroots(c8(a[ic]),c8(b[ic]),c8(c[ic]),byref(cd_r1),byref(cd_r2)) + r1[2] = cd_r2.value + r2[2] = cd_r1.value + + print(a[ic],b[ic],c[ic]) + print(r1) + print(r2) + + #PlotQuadAndRoots(a[ic],b[ic],c[ic],r1,r2) + + +def PlotQuadAndRoots(a,b,c,d,r1,r2): + + fig, axs = plt.subplots(ncols=1,nrows=1,figsize=(8,8)) + ax1s = axs.reshape(-1) + ic=0 + + npts = 1000 + + for i in range(npts): + print(i) + + + +# code.interact(local=dict(globals(), **locals())) + +# Helper code to plot negative logs + +def semilogneg(x): + + y = np.sign(x)*np.log(abs(x)) + return(y) + +def semilog10net(x): + + y = np.sign(x)*np.log10(abs(x)) + return(y) + + +# ======================================================================================= +# This is the actual call to main + +if __name__ == "__main__": + main(sys.argv) diff --git a/functional_unit_testing/math_utils/bld/README b/functional_unit_testing/math_utils/bld/README new file mode 100644 index 0000000000..4e67e5f091 --- /dev/null +++ b/functional_unit_testing/math_utils/bld/README @@ -0,0 +1 @@ +PLACEHOLDER FOR DIR \ No newline at end of file diff --git a/functional_unit_testing/math_utils/build_math_objects.sh b/functional_unit_testing/math_utils/build_math_objects.sh new file mode 100755 index 0000000000..40ac3eb9d1 --- /dev/null +++ b/functional_unit_testing/math_utils/build_math_objects.sh @@ -0,0 +1,47 @@ +#!/bin/bash + +# Path to FATES src + +FC='gfortran' + +F_OPTS="-shared -fPIC -g -ffpe-trap=zero,overflow,underflow -fbacktrace -fbounds-check" +#F_OPTS="-shared -fPIC -O" + + +MOD_FLAG="-J" + +rm -f bld/*.o +rm -f bld/*.mod + + +# First copy over the FatesConstants file, but change the types of the fates_r8 and fates_int + +old_fates_r8_str=`grep -e integer ../../main/FatesConstantsMod.F90 | grep fates_r8 | sed 's/^[ \t]*//;s/[ \t]*$//'` +new_fates_r8_str='use iso_c_binding, only: fates_r8 => c_double' + +old_fates_int_str=`grep -e integer ../../main/FatesConstantsMod.F90 | grep fates_int | sed 's/^[ \t]*//;s/[ \t]*$//'` +new_fates_int_str='use iso_c_binding, only: fates_int => c_int' + +# Add the new lines (need position change, don't swap) + +sed "/implicit none/i $new_fates_r8_str" ../../main/FatesConstantsMod.F90 > f90_src/FatesConstantsMod.F90 +sed -i "/implicit none/i $new_fates_int_str" f90_src/FatesConstantsMod.F90 +sed -i "/private /i public :: fates_r8" f90_src/FatesConstantsMod.F90 +sed -i "/private /i public :: fates_int" f90_src/FatesConstantsMod.F90 + +# Delete the old lines + +sed -i "/$old_fates_r8_str/d" f90_src/FatesConstantsMod.F90 +sed -i "/$old_fates_int_str/d" f90_src/FatesConstantsMod.F90 + +# Build the new file with constants + +${FC} ${F_OPTS} -I bld/ ${MOD_FLAG} bld/ -o bld/FatesConstantsMod.o f90_src/FatesConstantsMod.F90 + +${FC} ${F_OPTS} -I bld/ ${MOD_FLAG} bld/ -o bld/UnitWrapMod.o f90_src/UnitWrapMod.F90 + +${FC} ${F_OPTS} -I bld/ ${MOD_FLAG} bld/ -o bld/FatesUtilsMod.o ../../main/FatesUtilsMod.F90 + + + + diff --git a/functional_unit_testing/math_utils/f90_src/UnitWrapMod.F90 b/functional_unit_testing/math_utils/f90_src/UnitWrapMod.F90 new file mode 100644 index 0000000000..f12311655a --- /dev/null +++ b/functional_unit_testing/math_utils/f90_src/UnitWrapMod.F90 @@ -0,0 +1,49 @@ + +! ======================================================================================= +! +! This file is an alternative to key files in the fates +! filesystem. Noteably, we replace fates_r8 and fates_in +! with types that work with "ctypes". This is +! a key step in working with python +! +! We also wrap FatesGlobals to reduce the dependancy +! cascade that it pulls in from shr_log_mod. +! +! ======================================================================================= + +module shr_log_mod + + use iso_c_binding, only : c_char + use iso_c_binding, only : c_int + + contains + + function shr_log_errMsg(source, line) result(ans) + character(kind=c_char,len=*), intent(in) :: source + integer(c_int), intent(in) :: line + character(kind=c_char,len=128) :: ans + + ans = "source: " // trim(source) // " line: " + end function shr_log_errMsg + +end module shr_log_mod + + +module FatesGlobals + + contains + + integer function fates_log() + fates_log = -1 + end function fates_log + + subroutine fates_endrun(msg) + + implicit none + character(len=*), intent(in) :: msg ! string to be printed + + stop + + end subroutine fates_endrun + +end module FatesGlobals diff --git a/main/FatesUtilsMod.F90 b/main/FatesUtilsMod.F90 index 4699a6aa60..be955d5da0 100644 --- a/main/FatesUtilsMod.F90 +++ b/main/FatesUtilsMod.F90 @@ -5,7 +5,11 @@ module FatesUtilsMod use FatesConstantsMod, only : r8 => fates_r8 use FatesGlobals, only : fates_log - + use FatesConstantsMod, only : nearzero + use FatesGlobals, only : endrun => fates_endrun + + use shr_log_mod , only : errMsg => shr_log_errMsg + implicit none private ! Modules are private by default @@ -14,7 +18,12 @@ module FatesUtilsMod public :: check_var_real public :: GetNeighborDistance public :: FindIndex - + public :: QuadraticRootsNSWC + public :: QuadraticRootsSridharachary + + character(len=*), parameter, private :: sourcefile = & + __FILE__ + contains @@ -176,6 +185,97 @@ function FindIndex(input_string_array,string_to_match) result(array_index) end do end function FindIndex - + + subroutine QuadraticRootsNSWC(a,b,c,root1,root2) + + ! This code is based off of routines from the NSWC Mathematics Subroutine Library + ! From the NSWC README (https://github.com/jacobwilliams/nswc) + ! "The NSWC Mathematics Subroutine Library is a collection of Fortran 77 routines + ! specializing in numerical mathematics collected and developed by the U.S. + ! Naval Surface Warfare Center. This software is made available, without cost, + ! to the general scientific community." + ! The F77 code was updated to modern fortran by Jacob Williams: + ! https://jacobwilliams.github.io/polyroots-fortran + ! The FATES adaptation of this aborts if only imaginary roots are generated + + + real(r8),intent(in) :: a , b , c !! coefficients + real(r8),intent(out) :: root1 ! sr !! real part of first root + real(r8),intent(out) :: root2 ! lr !! real part of second root + + real(r8) :: b1, d, e + + if ( abs(a)nearzero ) then + ! compute discriminant avoiding overflow + b1 = b/2.0_r8 + if ( abs(b1)=0.0_r8 ) d = -d + root1 = (-b1+d)/a + root2 = 0.0_r8 + if ( root1/=0.0_r8 ) root2 = (c/root1)/a + endif + else + root2 = 0.0_r8 + root1 = -b/a + endif + + end subroutine QuadraticRootsNSWC + + subroutine QuadraticRootsSridharachary(a,b,c,root1,root2) + + + real(r8),intent(in) :: a , b , c !! coefficients + real(r8),intent(out) :: root1 ! sr !! real part of first root + real(r8),intent(out) :: root2 ! lr !! real part of second root + real(r8) :: d ! discriminant + real(r8) :: das ! sqrt(abs(d)) + + ! If a is 0, then equation is not quadratic, but linear + if (abs(a) < nearzero ) then + root2 = 0.0_r8 + if ( abs(b)>nearzero ) root2 = -c/b + root1 = 0.0_r8 + return + end if + + d = b * b - 4._r8 * a * c + das = sqrt(abs(d)) + + if (d > nearzero) then + + root1 = (-b + das) / (2._r8 * a) + root2 = (-b - das) / (2._r8 * a) + + elseif (abs(d) <= nearzero) then + + root1 = -b / (2._r8 * a) + root2 = root1 + else + + write (fates_log(),*)'error, imaginary roots detected in quadratic solve' + call endrun(msg=errMsg(sourcefile, __LINE__)) + + end if + + end subroutine QuadraticRootsSridharachary + ! ====================================================================================== end module FatesUtilsMod