Commit 9b2b6136 authored by Graham, Aaron's avatar Graham, Aaron
Browse files

Merge branch 'function_tabulation' into 'master' 

Add FunctionTableGenerator

See merge request futility/Futility!361
parents 2ec09ea4 a9709d69

src/CMakeLists.txt

+1 −0
Original line number Diff line number Diff line
@@ -103,6 +103,7 @@ TRIBITS_ADD_LIBRARY(Utils 
      SchemaParser.f90

      SpeciesElements.f90

      FutilityComputingEnvironment.f90

      FunctionTableGenerator.f90

      NonLinearSolver.f90

      SolutionAcceleration.f90

      AndersonAcceleration.f90

src/ExtendedMath.f90

+82 −2
Original line number Diff line number Diff line
@@ -22,6 +22,7 @@ PUBLIC :: GreatestCommonDivisor 
PUBLIC :: LeastCommonMultiple

PUBLIC :: ATAN2PI

PUBLIC :: RotateQtrClockwise

PUBLIC :: ThirdOrderBickley



INTERFACE LeastCommonMultiple

  MODULE PROCEDURE LeastCommonMultiple_scalar
@@ -279,4 +280,83 @@ ELEMENTAL SUBROUTINE RotateQtrClockwise(x,y,nrot) 
  ENDSELECT

ENDSUBROUTINE RotateQtrClockwise

!

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

!> @brief It calculates third order Bickley function

!> @param x The argument of the third order Bickley function

!> @returns y The result of the third order Bickley function

!>

PURE FUNCTION ThirdOrderBickley(x) RESULT(y)

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

  REAL(SRK) :: y



  !local variables

  REAL(SRK) :: xx



  xx=ABS(x)

  IF(xx < 0.05_SRK) THEN !Most common case

    y=(0.7266088_SRK*xx-0.9990226_SRK)*xx+0.7853961_SRK

  ELSEIF(xx > 9.0_SRK) THEN !Second most common case

    y=0.0_SRK

  ELSEIF(xx < 0.10_SRK) THEN

    y=(0.6466375_SRK*xx-0.9912340_SRK)*xx+0.7852024_SRK

  ELSEIF(xx < 0.15_SRK) THEN

    y=(0.5856605_SRK*xx-0.9791293_SRK)*xx+0.7845986_SRK

  ELSEIF(xx < 0.20_SRK) THEN

    y=(0.5346648_SRK*xx-0.9638914_SRK)*xx+0.7834577_SRK

  ELSEIF(xx < 0.25_SRK) THEN

    y=(0.4907827_SRK*xx-0.9463843_SRK)*xx+0.7817094_SRK

  ELSEIF(xx < 0.30_SRK) THEN

    y=(0.4521752_SRK*xx-0.9271152_SRK)*xx+0.7793031_SRK

  ELSEIF(xx < 0.35_SRK) THEN

    y=(0.4177388_SRK*xx-0.9064822_SRK)*xx+0.7762107_SRK

  ELSEIF(xx < 0.40_SRK) THEN

    y=(0.3869945_SRK*xx-0.8849865_SRK)*xx+0.7724519_SRK

  ELSEIF(xx < 0.45_SRK) THEN

    y=(0.3590753_SRK*xx-0.8626685_SRK)*xx+0.7679903_SRK

  ELSEIF(xx < 0.50_SRK) THEN

    y=(0.3338676_SRK*xx-0.8400133_SRK)*xx+0.7628988_SRK

  ELSEIF(xx < 0.60_SRK) THEN

    y=(0.2998569_SRK*xx-0.8054172_SRK)*xx+0.7540982_SRK

  ELSEIF(xx < 0.70_SRK) THEN

    y=(0.2609154_SRK*xx-0.7587821_SRK)*xx+0.7401279_SRK

  ELSEIF(xx < 0.80_SRK) THEN

    y=(0.2278226_SRK*xx-0.7125290_SRK)*xx+0.7239594_SRK

  ELSEIF(xx < 0.90_SRK) THEN

    y=(0.1994999_SRK*xx-0.6672761_SRK)*xx+0.7058777_SRK

  ELSEIF(xx < 1.00_SRK) THEN

    y=(0.1751248_SRK*xx-0.6234536_SRK)*xx+0.6861762_SRK

  ELSEIF(xx < 1.40_SRK) THEN

    y=((-0.05337485_SRK*xx+0.3203223_SRK)*xx-0.7538355_SRK)*xx+0.7247294_SRK

  ELSEIF(xx < 1.80_SRK) THEN

    y=((-0.03146833_SRK*xx+0.2295280_SRK)*xx-0.6279752_SRK)*xx+0.6663720_SRK

  ELSEIF(xx < 2.20_SRK) THEN

    y=((-0.01906198_SRK*xx+0.1631667_SRK)*xx-0.5094124_SRK)*xx+0.5956163_SRK

  ELSEIF(xx < 2.60_SRK) THEN

    y=((-0.01174752_SRK*xx+0.1152418_SRK)*xx-0.4046007_SRK)*xx+0.5191031_SRK

  ELSEIF(xx < 3.0_SRK) THEN

    y=((-0.007328415_SRK*xx+0.08097913_SRK)*xx-0.3159648_SRK)*xx+0.4425954_SRK

  ELSEIF(xx < 3.40_SRK) THEN

    y=((-0.004617254_SRK*xx+0.05669960_SRK)*xx-0.2434341_SRK)*xx+0.3703178_SRK

  ELSEIF(xx < 3.80_SRK) THEN

    y=(0.007923547_SRK*xx-0.07158569_SRK)*xx+0.1684022_SRK

  ELSEIF(xx < 4.20_SRK) THEN

    y=(0.005095111_SRK*xx-0.05016344_SRK)*xx+0.1278307_SRK

  ELSEIF(xx < 4.60_SRK) THEN

    y=(0.003286040_SRK*xx-0.03501524_SRK)*xx+0.09611422_SRK

  ELSEIF(xx < 5.00_SRK) THEN

    y=(0.002126242_SRK*xx-0.02437465_SRK)*xx+0.07170491_SRK

  ELSEIF(xx < 5.80_SRK) THEN

    y=(0.001123687_SRK*xx-0.01425519_SRK)*xx+0.04616317_SRK

  ELSEIF(xx < 6.60_SRK) THEN

    y=(4.762937e-4_SRK*xx-6.810124e-3_SRK)*xx+0.02475115_SRK

  ELSEIF(xx < 7.40_SRK) THEN

    y=(2.031843e-4_SRK*xx-3.232035e-3_SRK)*xx+0.01302864_SRK

  ELSEIF(xx < 8.20_SRK) THEN

    y=(8.701440e-5_SRK*xx-1.524126e-3_SRK)*xx+6.749972e-3_SRK

  ELSE !Range 8.2 <= xx <= 9.0, but the > 9.0 portion was already handled

    y=(3.742673e-5_SRK*xx-7.157367e-4_SRK)*xx+3.454768e-3_SRK

  ENDIF



ENDFUNCTION ThirdOrderBickley

!

ENDMODULE ExtendedMath

src/FunctionTableGenerator.f90

0 → 100644
+381 −0
Original line number Diff line number Diff line 
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!

!                          Futility Development Group                          !

!                             All rights reserved.                             !

!                                                                              !

! Futility is a jointly-maintained, open-source project between the University !

! of Michigan and Oak Ridge National Laboratory.  The copyright and license    !

! can be found in LICENSE.txt in the head directory of this repository.        !

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!

!> @brief This module defines the @c FunctionTable object and a factory subroutine

!> to generate the table

!>

!> The @c generateFunctionTable function should be called to retrieve a pointer

!> to the new @c FunctionTable.  The function to be tabluated is provided as an

!> argument.  A parameter list provides such parameters as the minimum and maximum

!> values for which the function should be tabulated, behavior for values outside

!> that range, and generation of the table.  The table can be generated using

!> fixed spacing, or it can be dynamically refined until client-specified criteria

!> such as maximum error are achieved.

!>

!> Once the @c FunctionTable pointer is returned to the client, that client is

!> now the owner of the pointer and the memory.  This module provides no

!> mechanisms for memory management of the table.  The client should call

!> @c evaluate on the table to obtain the value of the function for an input.

!> Once the client is done with the table, it should call @c clear, then

!> deallocate the pointer.

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!

MODULE FunctionTableGeneratorModule

#include "Futility_DBC.h"

USE Futility_DBC

USE IntrType

USE Constants_Conversion

USE IO_Strings

USE ParameterLists

USE FutilityComputingEnvironmentModule



IMPLICIT NONE

PRIVATE



PUBLIC :: FunctionTable

PUBLIC :: generateFunctionTable



#ifdef UNIT_TEST

PUBLIC :: generateTableData

PUBLIC :: runTestFunctionTableGenerator

INTERFACE

  MODULE SUBROUTINE runTestFunctionTableGenerator()

  ENDSUBROUTINE runTestFunctionTableGenerator

ENDINTERFACE

#endif



INTERFACE generateFunctionTable

  MODULE PROCEDURE generateFunctionTable1D

ENDINTERFACE generateFunctionTable



!> @brief abstract interface for a 1D tabulated function

!> @param x the input value

!> @returns y the output value

ABSTRACT INTERFACE

  FUNCTION tabulatedFunction1D(x) RESULT(y)

    IMPORT :: SRK

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

    REAL(SRK) :: y

  ENDFUNCTION tabulatedFunction1D

ENDINTERFACE



!> @brief Stores tabulated data for a function and returns the approximate function value

TYPE :: FunctionTable

  PRIVATE

  !> Logical indicating initialization status

  LOGICAL(SBK) :: isInit=.FALSE.

  !> Logical indicating if the table has been populated with data

  LOGICAL(SBK) :: hasData=.FALSE.

  !> Logical indicating if an error should be thrown for input values greater than @c maxRange

  LOGICAL(SBK) :: boundsErrorHigh=.TRUE.

  !> Logical indicating if an error should be thrown for input values less than @c minRange

  LOGICAL(SBK) :: boundsErrorLow=.TRUE.

  !> Number of data points contained in the table

  INTEGER(SIK) :: nPoints=0_SIK

  !> The minimum input value stored on the table

  REAL(SRK) :: minRange=HUGE(ZERO)

  !> The maximum input value stored on the table

  REAL(SRK) :: maxRange=-HUGE(ZERO)

  !> The return value for inputs greater than @c maxRange if @c boundsErrorHigh is false

  REAL(SRK) :: highDefaultValue=ZERO

  !> The return value for inputs greater than @c minRange if @c boundsErrorLow is false

  REAL(SRK) :: lowDefaultValue=ZERO

  !> The spacing between input values stored by the table

  REAL(SRK) :: spacing=-HUGE(ZERO)

  !> The inverse of @c spacing; stored for optimization purposes

  REAL(SRK) :: inverseSpacing=-HUGE(ZERO)

  !> The list of input values for which the function has been explicitly evaluated

  REAL(SRK),ALLOCATABLE :: funcCoordinates(:)

  !> The list of function values calculated from @c funcCoordinates and used for interpolation

  REAL(SRK),ALLOCATABLE :: funcValues(:)

  !> Pointer to the function used for tabulation

  PROCEDURE(tabulatedFunction1D),NOPASS,POINTER :: func => NULL()

  !> Pointer to the computing environment

  CLASS(FutilityComputingEnvironment),POINTER :: ce => NULL()

  !

  CONTAINS

    !> @copybrief FunctionTableGeneratorModule::initFunctionTable1D

    !> @copydetails FunctionTableGeneratorModule::initFunctionTable1D

    PROCEDURE,PASS :: init => initFunctionTable1D

    !> @copybrief FunctionTableGeneratorModule::clearFunctionTable1D

    !> @copydetails FunctionTableGeneratorModule::clearFunctionTable1D

    PROCEDURE,PASS :: clear => clearFunctionTable1D

    !> @copybrief FunctionTableGeneratorModule::evaluateFunctionTable1D

    !> @copydetails FunctionTableGeneratorModule::evaluateFunctionTable1D

    PROCEDURE,PASS :: evaluate => evaluateFunctionTable1D

ENDTYPE FunctionTable



CHARACTER(LEN=*),PARAMETER :: modName='FunctionTableGenerator'

!

!===============================================================================

CONTAINS

!

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

!> @brief generates a @c FunctionTable object and returns a pointer to it

!> @param ce the computing environment

!> @param params the input parameter list controlling the generation of the table

!> @param func the procedure to be tabulated; must have an interface matching @c tabulatedFunction1D

!> @returns table pointer to the initialized function table

!>

FUNCTION generateFunctionTable1D(ce,params,func) RESULT(table)

  CLASS(FutilityComputingEnvironment),TARGET,INTENT(INOUT) :: ce

  CLASS(ParamType),INTENT(IN) :: params

  PROCEDURE(tabulatedFunction1D) :: func

  CLASS(FunctionTable),POINTER :: table

  !

  REAL(SRK) :: maxRelError,maxRelErrorData,maxAbsError,maxAbsErrorData



  REQUIRE(params%has('spacing') /= (params%has('maxRelError') .OR. params%has('maxAbsError')))



  ALLOCATE(table)

  CALL table%init(ce,params,func)



  !Now setup the table.  If spacing is specified, just call setup once for that spacing

  IF(params%has('spacing')) THEN

    CALL params%get('spacing',table%spacing)

    REQUIRE(table%spacing > ZERO)

    CALL generateTableData(table)

  !If error is specified, start with 1000 points and refine until the error is small enough

  ELSE

    IF(params%has('maxRelError')) THEN

      CALL params%get('maxRelError',maxRelError)

      REQUIRE(maxRelError > ZERO)

    ELSE

      maxRelError=HUGE(maxRelError)

    ENDIF

    IF(params%has('maxAbsError')) THEN

      CALL params%get('maxAbsError',maxAbsError)

      REQUIRE(maxAbsError > ZERO)

    ELSE

      maxAbsError=HUGE(maxAbsError)

    ENDIF

    table%spacing = (table%maxRange - table%minRange)*0.001_SRK

    maxRelErrorData = HUGE(maxRelError)

    maxAbsErrorData = HUGE(maxAbsError)

    DO WHILE(maxRelErrorData > maxRelError .OR. maxAbsErrorData > maxAbsError)

      CALL generateTableData(table)

      CALL evaluateTableQuality(table,maxRelErrorData,maxAbsErrorData)

    ENDDO

  ENDIF



ENDFUNCTION generateFunctionTable1D

!

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

!> @brief Initializes a function table without populating its data

!> @param this the @c FunctionTable to initialize

!> @param ce the computing environment

!> @param params the input parameter list controlling the initialization of the table

!> @param func the procedure to be tabulated; must have an interface matching @c tabulatedFunction1D

!>

SUBROUTINE initFunctionTable1D(this,ce,params,func)

  CLASS(FunctionTable),INTENT(INOUT) :: this

  CLASS(FutilityComputingEnvironment),TARGET,INTENT(INOUT) :: ce

  CLASS(ParamType),INTENT(IN) :: params

  PROCEDURE(tabulatedFunction1D) :: func



  REQUIRE(params%has('minRange'))

  REQUIRE(params%has('maxRange'))

  REQUIRE(params%has('boundsErrorHigh'))

  REQUIRE(params%has('boundsErrorLow'))



  !Set the boundaries and default values (or error behavior) outside the boundaries

  CALL params%get('boundsErrorHigh',this%boundsErrorHigh)

  CALL params%get('boundsErrorLow',this%boundsErrorLow)

  IF(.NOT.this%boundsErrorHigh) THEN

    IF(params%has('highDefaultValue')) THEN

      CALL params%get('highDefaultValue',this%highDefaultValue)

    ENDIF

  ENDIF

  IF(.NOT.this%boundsErrorLow) THEN

    IF(params%has('lowDefaultValue')) THEN

      CALL params%get('lowDefaultValue',this%lowDefaultValue)

    ENDIF

  ENDIF

  CALL params%get('minRange',this%minRange)

  CALL params%get('maxRange',this%maxRange)

  REQUIRE(this%minRange < this%maxRange)



  !Associate pointers

  this%func => func

  this%ce => ce



  this%isInit=.TRUE.

  this%hasData=.FALSE.



ENDSUBROUTINE initFunctionTable1D

!

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

!> @brief Clears a @c FunctionTable object

!> @param this the table to clear

!>

SUBROUTINE clearFunctionTable1D(this)

  CLASS(FunctionTable),INTENT(INOUT) :: this



  this%isInit = .FALSE.

  this%hasData = .FALSE.

  this%boundsErrorHigh = .TRUE.

  this%boundsErrorLow = .TRUE.

  this%nPoints = 0_SIK

  this%minRange = HUGE(ZERO)

  this%maxRange = -HUGE(ZERO)

  this%highDefaultValue = ZERO

  this%lowDefaultValue = ZERO

  this%spacing = -HUGE(ZERO)

  this%inverseSpacing = -HUGE(ZERO)

  IF(ALLOCATED(this%funcCoordinates)) DEALLOCATE(this%funcCoordinates)

  IF(ALLOCATED(this%funcValues)) DEALLOCATE(this%funcValues)

  this%func => NULL()

  this%ce => NULL()



ENDSUBROUTINE clearFunctionTable1D

!

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

!> @brief Approximately evaluates the tabulated function at @c x

!> @param this the @c FunctionTable object

!> @param x the input value for the function

!> @returns y the approximate function values at @c x

!>

FUNCTION evaluateFunctionTable1D(this,x) RESULT(y)

  CLASS(FunctionTable),INTENT(IN) :: this

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

  REAL(SRK) :: y

  !

  CHARACTER(LEN=*),PARAMETER :: myName='evaluateFunctionTable1D'

  INTEGER :: tableIndex

  REAL(SRK) :: interval, fraction



  REQUIRE(this%isInit)

  REQUIRE(this%hasData)

  REQUIRE(this%inverseSpacing > ZERO)

  REQUIRE(ALLOCATED(this%funcCoordinates))

  REQUIRE(ALLOCATED(this%funcValues))



  interval = x - this%minRange

  tableIndex = FLOOR(interval * this%inverseSpacing) + 1

  IF(tableIndex < 1) THEN

    IF(this%boundsErrorLow) THEN

      CALL this%ce%exceptHandler%raiseError(modName//'::'//myName// &

          ' - Value '//str(x)//' is below the minimum table value '//str(this%minRange)//'!')

    ELSE

      y = this%lowDefaultValue

    ENDIF

  ELSEIF(tableIndex >= this%nPoints) THEN !Use >= to avoid index out of bounds in ELSE statement

    IF(this%boundsErrorHigh) THEN

      CALL this%ce%exceptHandler%raiseError(modName//'::'//myName// &

          ' - Value '//str(x)//' is above the maximum table value '//str(this%maxRange)//'!')

    ELSE

      y = this%highDefaultValue

    ENDIF

  ELSE

    fraction = (x - this%funcCoordinates(tableIndex)) * this%inverseSpacing

    y = (ONE - fraction) * this%funcValues(tableIndex) + fraction * this%funcValues(tableIndex+1)

  ENDIF



ENDFUNCTION evaluateFunctionTable1D

!

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

!> @brief Generates the data for a @c FunctionTable object

!> @param table the object to populate data for

!>

!> If @c hasData is true, then it is assumed that the function table data is being

!> refined.  The spacing will be cut in half and the data regenerated.  If false, then

!> new data will be populated using whatever the current value of @c spacing is.

!>

SUBROUTINE generateTableData(table)

  CLASS(FunctionTable),INTENT(INOUT) :: table

  !

  INTEGER(SIK) :: i,stride

  REAL(SRK) :: x

  REAL(SRK),ALLOCATABLE :: oldValues(:),oldCoordinates(:)



  REQUIRE(ASSOCIATED(table%func))

  REQUIRE(table%maxRange > -HUGE(table%maxRange))

  REQUIRE(table%minRange < HUGE(table%minRange))

  REQUIRE(table%spacing > ZERO)



  IF(table%hasData) THEN

    !Assume we've already generated the table once and are refining by 2x

    IF((SIZE(table%funcValues) > 2) .AND. (MOD(SIZE(table%funcValues),2) == 1)) THEN

      CALL MOVE_ALLOC(table%funcValues,oldValues)

      CALL MOVE_ALLOC(table%funcCoordinates,oldCoordinates)

      stride = 2

      table%spacing = table%spacing * HALF

    !Generated some table that doesn't conform to how we're going to generate it this time,

    ! so just throw out the old values

    ELSE

      DEALLOCATE(table%funcValues)

      DEALLOCATE(table%funcCoordinates)

      stride = 1

    ENDIF

  ELSE

    stride = 1

  ENDIF



  table%inverseSpacing = ONE / table%spacing

  table%nPoints = NINT((table%maxRange - table%minRange) * table%inverseSpacing) + 1

  ALLOCATE(table%funcValues(table%nPoints))

  ALLOCATE(table%funcCoordinates(table%nPoints))





  IF(stride == 2) THEN

    table%funcValues(1:table%nPoints:2)=oldValues(:)

    table%funcCoordinates(1:table%nPoints:2)=oldCoordinates(:)

    DEALLOCATE(oldValues)

    DEALLOCATE(oldCoordinates)

  ENDIF



  DO i=stride,table%nPoints,stride

    x = table%minRange + REAL(i-1,SRK)*table%spacing

    table%funcCoordinates(i) = x

    table%funcValues(i) = table%func(x)

  ENDDO !i



  table%hasData=.TRUE.



ENDSUBROUTINE generateTableData

!

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

!> @brief Evaluates the quality of a @c FunctionTable compared to the real function

!> @param table the @c FunctionTable to evaluate

!> @returns maxRelError the maximum relative error of the function table

!> @returns maxAbsError the maximum absolute error of the function table

!>

!> If the function returns 0.0 within double precision tolerance, then the value

!> of the table at that input is used to calculate the relative error.  This may

!> be of importance when generating the table dynamically.

!>

SUBROUTINE evaluateTableQuality(table,maxRelError,maxAbsError)

  CLASS(FunctionTable),INTENT(IN) :: table

  REAL(SRK),INTENT(OUT) :: maxRelError

  REAL(SRK),INTENT(OUT) :: maxAbsError

  !

  INTEGER(SIK) :: i

  REAL(SRK) :: x,yTab,yRef



  REQUIRE(ALLOCATED(table%funcValues))

  REQUIRE(ASSOCIATED(table%func))

  REQUIRE(table%maxRange > -HUGE(table%maxRange))

  REQUIRE(table%minRange < HUGE(table%minRange))

  REQUIRE(table%spacing > ZERO)



  maxRelError = ZERO

  maxAbsError = ZERO

  DO i=1,table%nPoints-1

    x = table%minRange + (REAL(i-1,SRK)+HALF) * table%spacing

    yTab = table%evaluate(x)

    yRef = table%func(x)

    IF((yRef .APPROXEQA. ZERO) .AND. (.NOT.(yTab .APPROXEQA. ZERO))) THEN

      maxRelError = MAX(maxRelError,yTab)

    ELSE

      maxRelError = MAX(maxRelError,ABS((yTab - yRef)/yRef))

    ENDIF

    maxAbsError = MAX(maxAbsError,ABS(yTab - yRef))

  ENDDO



ENDSUBROUTINE evaluateTableQuality

!

ENDMODULE FunctionTableGeneratorModule
 
 No newline at end of file

unit_tests/CMakeLists.txt

+1 −0
Original line number Diff line number Diff line
@@ -76,6 +76,7 @@ SET(UNIT_TEST_NAMES 
    testSchemaParser

    testSpeciesElements

    testRSORprecon

    testFunctionTableGenerator

)



IF("${PROJECT_NAME}" STREQUAL "Futility" OR "${PROJECT_NAME}" STREQUAL "MPACT")

unit_tests/testExtendedMath/testExtendedMath.f90

+69 −0

File changed.

Preview size limit exceeded, changes collapsed.