Commit a0c70634 authored by aarograh's avatar aarograh Committed by Henderson, Shane
Browse files

Implement the Regularized Incomplete Beta Function 

Squash branch 'beta_cdf' into 'master'

* Fix use of uninit variable in REQUIRE

* Use SOFTEQ for beta function test

* Fix unused variables

* Refine the naming for the beta function

* Fix missing FINALIZE_TEST and precision failure in testWaterSatProperties

* Implement betaCDF function & tests

This is being used in some new rod cusping treatments

**Developer Checklist:**
- [x] Have you done a self-review after creating the merge request?
- [x] Have you filled in the Merge Request information (title, description) thoroughly?
- [x] Have you updated the relevant tickets (if this MR is linked to any VERA-dev tickets)?
- [x] Have you addressed all suggested feedback and commented on it to let the reviewer know? (Do not resolve discussions that the reviewer started)

**Reviewer Checklist:**
- [x] Have you confirmed all discussions were adequately addressed and resolved them all?
- [x] Does it conform to formatting guidelines?
- [x] Are there adequate and clear comments?
- [x] Is the design clean and sensible?
- [x] Are the changes optimal/efficient?
- [x] Were sufficient DBC checks added?
- [x] Are there unit tests? (if necessary)
- [x] Is the MR description clear, including a link to the VERA-Dev issue if appropriate?

**PSM Checklist**
- [x] Have you confirmed that all discussions were addressed, or that follow-on issues have been created for them?
- [x] Have you confirmed sufficient testing was conducted?
- [x] Does this impact other repositories?
- [x] Does the MR have an adequate description?
- [x] If the MR has multiple commits, did you set the MR to squash merge?

See merge request https://code.ornl.gov/futility/Futility/-/merge_requests/409
parent b846ca04

src/ExtendedMath.f90

+70 −0
Original line number Diff line number Diff line
@@ -11,6 +11,8 @@ 
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!

MODULE ExtendedMath



#include "Futility_DBC.h"

USE Futility_DBC

USE IntrType

USE Constants_Conversion

IMPLICIT NONE
@@ -23,6 +25,7 @@ PUBLIC :: LeastCommonMultiple 
PUBLIC :: ATAN2PI

PUBLIC :: RotateQtrClockwise

PUBLIC :: ThirdOrderBickley

PUBLIC :: RegularizedIncompleteBeta



INTERFACE LeastCommonMultiple

  MODULE PROCEDURE LeastCommonMultiple_scalar
@@ -360,4 +363,71 @@ PURE FUNCTION ThirdOrderBickley(x) RESULT(y) 


ENDFUNCTION ThirdOrderBickley

!

!-------------------------------------------------------------------------------

!> @brief Computes the regularized incomplete beta function

!> @param val the value at which to evaluate the function

!> @param a the first shape parameter of the function

!> @param b the second shape parameter of the function

!> @param width the width of the integration step (optional, defaults to 0.001)

!> @returns y the value of the function at val

!>

!> Ranges for the a, b parameters and the resulting shape:

!>   a = b = 1.0: uniform distribution

!>   a > 1, b > 1: Bell-shaped (symmetric if a = b)

!>   a < 1, b < 1: U-shaped (infinite density at 0 and 1)

!>   a < 1, b > 1: Skewed toward 0

!>   a > 1, b < 1: Skewed toward 1

!>   a >> b: Very steep rise near 0, flat near 1

!>   a << b: Very flat near 0, steep rise near 1

!>

FUNCTION RegularizedIncompleteBeta(val,a,b,width) RESULT(y)

  REAL(SRK),INTENT(IN) :: val

  REAL(SRK),INTENT(IN) :: a

  REAL(SRK),INTENT(IN) :: b

  REAL(SRK),INTENT(IN),OPTIONAL :: width

  REAL(SRK) :: y

  !

  INTEGER(SIK) :: i,intervals

  REAL(SRK) :: x,dx,x1,x2,y1,y2



  REQUIRE(val >= 0.0_SRK)

  REQUIRE(val <= 1.0_SRK)

  REQUIRE(a > 0.0_SRK)

  REQUIRE(b > 0.0_SRK)

  dx = 0.001_SRK

  IF(PRESENT(width)) THEN

    dx = width

  ENDIF

  REQUIRE(dx > 0.0_SRK)

  REQUIRE(dx <= 1.0_SRK)



  ! Trouble evaluation at the end points, but the function goes to 0 at 0 and 1 at 1,

  ! so just early return if we're close enough to the end points

  x = val

  IF((x .APPROXEQ. ZERO) .OR. (x .APPROXEQ. ONE)) THEN

    y = x

    RETURN

  ENDIF



  ! Numerically integrate the function from 0 to x

  y = 0.0_SRK

  x1 = 1.0E-14_SRK

  y1 = 0.0_SRK

  intervals = CEILING((x - x1) / dx)

  dx = (x - x1) / intervals

  DO i = 1, intervals

    x2 = x1 + dx

    y2 = x2**(a - 1.0_SRK) * (1.0_SRK - x2)**(b - 1.0_SRK)

    y = y + 0.5_SRK * (y1 + y2) * dx

    x1 = x2

    y1 = y2

  ENDDO !i



  ! Normalize the result by the "regularized" beta function, which is just

  ! BetaCDF(1.0_SRK,a,b), which is also equal to

  ! GAMMF(a) * GAMMAF(b) / GAMMAF(a + b)

  y = y * GAMMAF(a + b) / (GAMMAF(a) * GAMMAF(b))



END FUNCTION RegularizedIncompleteBeta

!

ENDMODULE ExtendedMath

unit_tests/testExtendedMath/testExtendedMath.f90

+10 −0
Original line number Diff line number Diff line
@@ -24,6 +24,7 @@ REGISTER_SUBTEST('GCD',testLCM) 
REGISTER_SUBTEST('ATAN2PI',testATAN2PI)

REGISTER_SUBTEST('RotateQtrClockwise',testRotateQtrClockwise)

REGISTER_SUBTEST("testThirdOrderBickley", testThirdOrderBickley)

REGISTER_SUBTEST("testRegularizedIncompleteBeta", testRegularizedIncompleteBeta)



FINALIZE_TEST()

!
@@ -316,4 +317,13 @@ SUBROUTINE testThirdOrderBickley() 
!8.600     0.000067544575386678789251476164807145393151299089940651099129594815501432040017895

ENDSUBROUTINE testThirdOrderBickley

!

!------------------------------------------------------------------------------

SUBROUTINE testRegularizedIncompleteBeta()



  ASSERT_APPROXEQ(regularizedincompletebeta(0.0_SRK, 2.53_SRK, 0.15_SRK, 0.000001_SRK), 0.0_SRK, 'regularizedincompletebeta(0.0)')

  ASSERT_APPROXEQ(regularizedincompletebeta(1.0_SRK, 2.53_SRK, 0.15_SRK, 0.000001_SRK), 1.0_SRK, 'regularizedincompletebeta(1.0)')

  ASSERT_APPROXEQ(regularizedincompletebeta(0.99_SRK, 2.53_SRK, 0.15_SRK, 0.000001_SRK), 0.39979661542029937_SRK, 'regularizedincompletebeta(0.99)')



ENDSUBROUTINE testRegularizedIncompleteBeta

!

ENDPROGRAM testExtendedMath

unit_tests/testWaterSatProperties/testWaterSatProperties.f90

+32 −30
Original line number Diff line number Diff line
@@ -21,56 +21,58 @@ CREATE_TEST('WATERSATPROPERTIES') 
!Test uninit

COMPONENT_TEST('Uninititialized')

refval=-HUGE(refval)

ASSERT(WaterSatProperties_GetTemp(1040.0_SRK) == refval,'GetTemp')

ASSERT(WaterSatProperties_GetPres(1040.0_SRK) == refval,'GetPres')

ASSERT(WaterSatProperties_GetLiqDens(P=1040.0_SRK) == refval,'GetLiqDens(P)')

ASSERT(WaterSatProperties_GetLiqDens(T=560.0_SRK) == refval,'GetLiqDens(T)')

ASSERT(WaterSatProperties_GetVapDens(P=1040.0_SRK) == refval,'GetVapDens(P)')

ASSERT(WaterSatProperties_GetVapDens(T=560.0_SRK) == refval,'GetVapDens(T)')

ASSERT_EQ(WaterSatProperties_GetTemp(1040.0_SRK) , refval, 'GetTemp')

ASSERT_EQ(WaterSatProperties_GetPres(1040.0_SRK) , refval, 'GetPres')

ASSERT_EQ(WaterSatProperties_GetLiqDens(P=1040.0_SRK) , refval, 'GetLiqDens(P)')

ASSERT_EQ(WaterSatProperties_GetLiqDens(T=560.0_SRK) , refval, 'GetLiqDens(T)')

ASSERT_EQ(WaterSatProperties_GetVapDens(P=1040.0_SRK) , refval, 'GetVapDens(P)')

ASSERT_EQ(WaterSatProperties_GetVapDens(T=560.0_SRK) , refval, 'GetVapDens(T)')



CALL WaterSatProperties_Init()



!Nominal BWR Conditions

COMPONENT_TEST('Nominal BWR Properties')

refval=560.610988281243_SRK

ASSERT(refval .APPROXEQ. WaterSatProperties_GetTemp(1040.0_SRK),'GetTemp')

ASSERT_SOFTEQ(refval , WaterSatProperties_GetTemp(1040.0_SRK), 1.0E-11_SRK, 'GetTemp')

refval=1040.0_SRK

ASSERT(refval .APPROXEQ. WaterSatProperties_GetPres(560.610988281243_SRK),'GetTemp')

ASSERT_SOFTEQ(refval , WaterSatProperties_GetPres(560.610988281243_SRK), 1.0E-11_SRK, 'GetTemp')

refval=0.736690682364633_SRK

ASSERT(refval .APPROXEQ. WaterSatProperties_GetLiqDens(P=1040.0_SRK),'GetLiqDens(P)')

ASSERT(refval .APPROXEQ. WaterSatProperties_GetLiqDens(T=560.610988281243_SRK),'GetLiqDens(T)')

ASSERT_APPROXEQ(refval , WaterSatProperties_GetLiqDens(P=1040.0_SRK), 'GetLiqDens(P)')

ASSERT_APPROXEQ(refval , WaterSatProperties_GetLiqDens(T=560.610988281243_SRK), 'GetLiqDens(T)')

refval=3.752512913761709E-02_SRK

ASSERT(refval .APPROXEQ. WaterSatProperties_GetVapDens(P=1040.0_SRK),'GetLiqDens(P)')

ASSERT_APPROXEQ(refval , WaterSatProperties_GetVapDens(P=1040.0_SRK), 'GetLiqDens(P)')

testval=WaterSatProperties_GetVapDens(T=WaterSatProperties_GetTemp(1040.0_SRK))

ASSERT(refval .APPROXEQ. testval,'GetLiqDens(T)')

ASSERT_APPROXEQ(refval , testval, 'GetLiqDens(T)')



!Test exact

COMPONENT_TEST('Exact Evaluation')

refval=300.0_SRK

ASSERT(refval == WaterSatProperties_GetTemp(5.12970461890E-01_SRK),'GetTemp exact')

ASSERT_EQ(refval , WaterSatProperties_GetTemp(5.12970461890E-01_SRK), 'GetTemp exact')

refval=5.12970461890E-01_SRK

ASSERT(refval == WaterSatProperties_GetPres(300.0_SRK),'GetPres exact')

ASSERT_EQ(refval , WaterSatProperties_GetPres(300.0_SRK), 'GetPres exact')

refval=9.96513027530E-01_SRK

ASSERT(refval .APPROXEQ. WaterSatProperties_GetLiqDens(P=5.12970461890E-01_SRK),'GetLiqDens(P)')

ASSERT(refval .APPROXEQ. WaterSatProperties_GetLiqDens(T=300.0_SRK),'GetLiqDens(T)')

ASSERT_APPROXEQ(refval , WaterSatProperties_GetLiqDens(P=5.12970461890E-01_SRK), 'GetLiqDens(P)')

ASSERT_APPROXEQ(refval , WaterSatProperties_GetLiqDens(T=300.0_SRK), 'GetLiqDens(T)')

refval=2.55896736840E-05_SRK

ASSERT(refval .APPROXEQ. WaterSatProperties_GetVapDens(P=5.12970461890E-01_SRK),'GetLiqDens(P)')

ASSERT(refval .APPROXEQ. WaterSatProperties_GetVapDens(T=300.0_SRK),'GetLiqDens(T)')

ASSERT_APPROXEQ(refval , WaterSatProperties_GetVapDens(P=5.12970461890E-01_SRK), 'GetLiqDens(P)')

ASSERT_APPROXEQ(refval , WaterSatProperties_GetVapDens(T=300.0_SRK), 'GetLiqDens(T)')



!Test out of bounds

COMPONENT_TEST('Out of Bounds')

refval=-HUGE(refval)

ASSERT(WaterSatProperties_GetTemp(0.01_SRK) == refval,'GetTemp below')

ASSERT(WaterSatProperties_GetTemp(3300.11_SRK) == refval,'GetTemp above')

ASSERT(WaterSatProperties_GetPres(273.0_SRK) == refval,'GetPres below')

ASSERT(WaterSatProperties_GetPres(678.0_SRK) == refval,'GetPres above')

ASSERT(WaterSatProperties_GetLiqDens(P=0.01_SRK) == refval,'GetLiqDens P below')

ASSERT(WaterSatProperties_GetLiqDens(P=3300.11_SRK) == refval,'GetLiqDens P above')

ASSERT(WaterSatProperties_GetLiqDens(T=273.0_SRK) == refval,'GetLiqDens T below')

ASSERT(WaterSatProperties_GetLiqDens(T=678.0_SRK) == refval,'GetLiqDens T above')

ASSERT(WaterSatProperties_GetVapDens(P=0.01_SRK) == refval,'GetVapDens P below')

ASSERT(WaterSatProperties_GetVapDens(P=3300.11_SRK) == refval,'GetVapDens P above')

ASSERT(WaterSatProperties_GetVapDens(T=273.0_SRK) == refval,'GetVapDens T below')

ASSERT(WaterSatProperties_GetVapDens(T=678.0_SRK) == refval,'GetVapDens T above')

ASSERT_EQ(WaterSatProperties_GetTemp(0.01_SRK) , refval, 'GetTemp below')

ASSERT_EQ(WaterSatProperties_GetTemp(3300.11_SRK) , refval, 'GetTemp above')

ASSERT_EQ(WaterSatProperties_GetPres(273.0_SRK) , refval, 'GetPres below')

ASSERT_EQ(WaterSatProperties_GetPres(678.0_SRK) , refval, 'GetPres above')

ASSERT_EQ(WaterSatProperties_GetLiqDens(P=0.01_SRK) , refval, 'GetLiqDens P below')

ASSERT_EQ(WaterSatProperties_GetLiqDens(P=3300.11_SRK) , refval, 'GetLiqDens P above')

ASSERT_EQ(WaterSatProperties_GetLiqDens(T=273.0_SRK) , refval, 'GetLiqDens T below')

ASSERT_EQ(WaterSatProperties_GetLiqDens(T=678.0_SRK) , refval, 'GetLiqDens T above')

ASSERT_EQ(WaterSatProperties_GetVapDens(P=0.01_SRK) , refval, 'GetVapDens P below')

ASSERT_EQ(WaterSatProperties_GetVapDens(P=3300.11_SRK) , refval, 'GetVapDens P above')

ASSERT_EQ(WaterSatProperties_GetVapDens(T=273.0_SRK) , refval, 'GetVapDens T below')

ASSERT_EQ(WaterSatProperties_GetVapDens(T=678.0_SRK) , refval, 'GetVapDens T above')



FINALIZE_TEST()

!

ENDPROGRAM testWaterSatProperties