DAMASK_EICMD/src/polynomials.f90

197 lines
6.4 KiB
Fortran
Raw Normal View History

2022-01-31 19:35:15 +05:30
!--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, KU Leuven
2022-12-07 22:59:03 +05:30
!> @brief Polynomial representation for variable data.
2022-01-31 19:35:15 +05:30
!--------------------------------------------------------------------------------------------------
module polynomials
use prec
use IO
use YAML_parse
use YAML_types
implicit none(type,external)
2022-01-31 19:35:15 +05:30
private
type, public :: tPolynomial
real(pReal), dimension(:), allocatable :: coef
real(pReal) :: x_ref = huge(0.0_pReal)
2022-01-31 19:35:15 +05:30
contains
procedure, public :: at => eval
end type tPolynomial
interface polynomial
module procedure polynomial_from_coef
2022-12-07 22:59:03 +05:30
module procedure polynomial_from_dict
2022-01-31 19:35:15 +05:30
end interface polynomial
public :: &
polynomial, &
polynomials_init
contains
!--------------------------------------------------------------------------------------------------
!> @brief Run self-test.
!--------------------------------------------------------------------------------------------------
subroutine polynomials_init()
print'(/,1x,a)', '<<<+- polynomials init -+>>>'; flush(IO_STDOUT)
call selfTest()
end subroutine polynomials_init
!--------------------------------------------------------------------------------------------------
2022-12-07 22:59:03 +05:30
!> @brief Initialize a polynomial from coefficients.
2022-01-31 19:35:15 +05:30
!--------------------------------------------------------------------------------------------------
2022-05-08 01:09:43 +05:30
pure function polynomial_from_coef(coef,x_ref) result(p)
2022-01-31 19:35:15 +05:30
2022-05-08 01:09:43 +05:30
real(pReal), dimension(0:), intent(in) :: coef
2022-01-31 19:35:15 +05:30
real(pReal), intent(in) :: x_ref
type(tPolynomial) :: p
2022-05-08 01:09:43 +05:30
p%coef = coef
2022-01-31 19:35:15 +05:30
p%x_ref = x_ref
end function polynomial_from_coef
!--------------------------------------------------------------------------------------------------
2022-12-07 22:59:03 +05:30
!> @brief Initialize a polynomial from a dictionary with coefficients.
2022-01-31 19:35:15 +05:30
!--------------------------------------------------------------------------------------------------
function polynomial_from_dict(dict,y,x) result(p)
type(tDict), intent(in) :: dict
character(len=*), intent(in) :: y, x
type(tPolynomial) :: p
real(pReal), dimension(:), allocatable :: coef
real(pReal) :: x_ref
2022-05-08 02:23:31 +05:30
integer :: i, o
character(len=1) :: o_s
2022-01-31 19:35:15 +05:30
allocate(coef(1),source=dict%get_asReal(y))
2022-01-31 19:35:15 +05:30
if (dict%contains(y//','//x)) then
x_ref = dict%get_asReal(x//'_ref')
coef = [coef,dict%get_asReal(y//','//x)]
end if
2022-05-08 02:34:06 +05:30
do o = 2,4
2022-05-08 02:23:31 +05:30
write(o_s,'(I0.0)') o
if (dict%contains(y//','//x//'^'//o_s)) then
x_ref = dict%get_asReal(x//'_ref')
coef = [coef,[(0.0_pReal,i=size(coef),o-1)],dict%get_asReal(y//','//x//'^'//o_s)]
2022-05-08 02:23:31 +05:30
end if
end do
2022-01-31 19:35:15 +05:30
p = Polynomial(coef,x_ref)
end function polynomial_from_dict
!--------------------------------------------------------------------------------------------------
2022-12-07 22:59:03 +05:30
!> @brief Evaluate a polynomial.
2022-05-08 23:27:25 +05:30
!> @details https://nvlpubs.nist.gov/nistpubs/jres/71b/jresv71bn1p11_a1b.pdf (eq. 1.2)
2022-01-31 19:35:15 +05:30
!--------------------------------------------------------------------------------------------------
pure function eval(self,x) result(y)
class(tPolynomial), intent(in) :: self
real(pReal), intent(in) :: x
real(pReal) :: y
2022-05-08 23:27:25 +05:30
integer :: o
2022-01-31 19:35:15 +05:30
y = self%coef(ubound(self%coef,1))
2022-05-08 23:27:25 +05:30
do o = ubound(self%coef,1)-1, 0, -1
2022-05-08 21:17:20 +05:30
#ifndef __INTEL_LLVM_COMPILER
2022-05-08 23:27:25 +05:30
y = y*(x-self%x_ref) +self%coef(o)
#else
2022-05-08 23:27:25 +05:30
y = IEEE_FMA(y,x-self%x_ref,self%coef(o))
#endif
2022-06-09 02:36:01 +05:30
end do
2022-01-31 19:35:15 +05:30
end function eval
!--------------------------------------------------------------------------------------------------
!> @brief Check correctness of polynomical functionality.
!--------------------------------------------------------------------------------------------------
subroutine selfTest()
2022-01-31 19:35:15 +05:30
type(tPolynomial) :: p1, p2
2022-05-08 02:34:06 +05:30
real(pReal), dimension(5) :: coef
integer :: i
2022-05-08 15:42:42 +05:30
real(pReal) :: x_ref, x, y
type(tDict), pointer :: dict
character(len=pStringLen), dimension(size(coef)) :: coef_s
2022-01-31 19:35:15 +05:30
character(len=pStringLen) :: x_ref_s, x_s, YAML_s
2022-01-31 19:35:15 +05:30
call random_number(coef)
call random_number(x_ref)
call random_number(x)
coef = coef*10_pReal -0.5_pReal
x_ref = x_ref*10_pReal -0.5_pReal
x = x*10_pReal -0.5_pReal
2022-05-08 15:42:42 +05:30
p1 = polynomial([coef(1)],x_ref)
if (dNeq(p1%at(x),coef(1))) error stop 'polynomial: eval(constant)'
2022-01-31 19:35:15 +05:30
p1 = polynomial(coef,x_ref)
if (dNeq(p1%at(x_ref),coef(1))) error stop 'polynomial: @ref'
do i = 1, size(coef_s)
write(coef_s(i),*) coef(i)
end do
2022-01-31 19:35:15 +05:30
write(x_ref_s,*) x_ref
write(x_s,*) x
YAML_s = 'C: '//trim(adjustl(coef_s(1)))//IO_EOL//&
'C,T: '//trim(adjustl(coef_s(2)))//IO_EOL//&
'C,T^2: '//trim(adjustl(coef_s(3)))//IO_EOL//&
2022-05-08 02:34:06 +05:30
'C,T^3: '//trim(adjustl(coef_s(4)))//IO_EOL//&
'C,T^4: '//trim(adjustl(coef_s(5)))//IO_EOL//&
2022-01-31 19:35:15 +05:30
'T_ref: '//trim(adjustl(x_ref_s))//IO_EOL
dict => YAML_parse_str_asDict(trim(YAML_s))
p2 = polynomial(dict,'C','T')
if (dNeq(p1%at(x),p2%at(x),1.0e-6_pReal)) error stop 'polynomials: init'
2022-05-08 15:42:42 +05:30
y = coef(1)+coef(2)*(x-x_ref)+coef(3)*(x-x_ref)**2+coef(4)*(x-x_ref)**3+coef(5)*(x-x_ref)**4
if (dNeq(p1%at(x),y,1.0e-6_pReal)) error stop 'polynomials: eval(full)'
2022-01-31 19:35:15 +05:30
YAML_s = 'C: 0.0'//IO_EOL//&
'C,T: '//trim(adjustl(coef_s(2)))//IO_EOL//&
'T_ref: '//trim(adjustl(x_ref_s))//IO_EOL
dict => YAML_parse_str_asDict(trim(YAML_s))
p1 = polynomial(dict,'C','T')
2022-05-08 02:34:06 +05:30
if (dNeq(p1%at(x_ref+x),-p1%at(x_ref-x),1.0e-10_pReal)) error stop 'polynomials: eval(linear)'
2022-01-31 19:35:15 +05:30
YAML_s = 'C: 0.0'//IO_EOL//&
'C,T^2: '//trim(adjustl(coef_s(3)))//IO_EOL//&
'T_ref: '//trim(adjustl(x_ref_s))//IO_EOL
dict => YAML_parse_str_asDict(trim(YAML_s))
p1 = polynomial(dict,'C','T')
2022-05-08 02:34:06 +05:30
if (dNeq(p1%at(x_ref+x),p1%at(x_ref-x),1e-10_pReal)) error stop 'polynomials: eval(quadratic)'
2022-05-08 15:42:42 +05:30
YAML_s = 'Y: '//trim(adjustl(coef_s(1)))//IO_EOL//&
'Y,X^3: '//trim(adjustl(coef_s(2)))//IO_EOL//&
2022-05-08 02:34:06 +05:30
'X_ref: '//trim(adjustl(x_ref_s))//IO_EOL
dict => YAML_parse_str_asDict(trim(YAML_s))
p1 = polynomial(dict,'Y','X')
if (dNeq(p1%at(x_ref+x)-coef(1),-(p1%at(x_ref-x)-coef(1)),1.0e-8_pReal)) error stop 'polynomials: eval(cubic)'
2022-05-08 02:34:06 +05:30
2022-05-08 15:42:42 +05:30
YAML_s = 'Y: '//trim(adjustl(coef_s(1)))//IO_EOL//&
'Y,X^4: '//trim(adjustl(coef_s(2)))//IO_EOL//&
2022-05-08 02:34:06 +05:30
'X_ref: '//trim(adjustl(x_ref_s))//IO_EOL
dict => YAML_parse_str_asDict(trim(YAML_s))
p1 = polynomial(dict,'Y','X')
if (dNeq(p1%at(x_ref+x),p1%at(x_ref-x),1.0e-6_pReal)) error stop 'polynomials: eval(quartic)'
2022-05-08 02:34:06 +05:30
2022-01-31 19:35:15 +05:30
end subroutine selfTest
end module polynomials