DAMASK_EICMD/src/polynomials.f90

196 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, &
polynomials_selfTest
2022-01-31 19:35:15 +05:30
contains
!--------------------------------------------------------------------------------------------------
!> @brief Run self-test.
!--------------------------------------------------------------------------------------------------
subroutine polynomials_init()
print'(/,1x,a)', '<<<+- polynomials init -+>>>'; flush(IO_STDOUT)
call polynomials_selfTest()
2022-01-31 19:35:15 +05:30
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
real(pREAL), dimension(0:), intent(in) :: coef
real(pREAL), intent(in) :: x_ref
2022-01-31 19:35:15 +05:30
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
2022-05-08 02:23:31 +05:30
integer :: i, o
character :: 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)) coef = [coef,dict%get_asReal(y//','//x)]
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)) &
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 do
2022-01-31 19:35:15 +05:30
if (size(coef) > 1) then
p = polynomial(coef,dict%get_asReal(x//'_ref'))
else
p = polynomial(coef,-huge(1.0_pREAL))
end if
2022-01-31 19:35:15 +05:30
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-01-31 19:35:15 +05:30
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))
do o = ubound(self%coef,1)-1, 0, -1
y = y*(x-self%x_ref) + self%coef(o)
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 polynomials_selfTest()
2022-01-31 19:35:15 +05:30
type(tPolynomial) :: p1, p2
real(pREAL), dimension(5) :: coef
integer :: i
real(pREAL) :: x_ref, x, y
type(tDict), pointer :: dict
2023-06-04 10:47:38 +05:30
character(len=pSTRLEN), dimension(size(coef)) :: coef_s
character(len=pSTRLEN) :: x_ref_s, x_s, YAML_s
2022-01-31 19:35:15 +05:30
2022-01-31 19:35:15 +05:30
call random_number(coef)
call random_number(x_ref)
call random_number(x)
coef = 10_pREAL*(coef-0.5_pREAL)
x_ref = 10_pREAL*(x_ref-0.5_pREAL)
x = 10_pREAL*(x-0.5_pREAL)
2022-01-31 19:35:15 +05:30
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'
y = coef(1)*(x-x_ref)**0 &
+ coef(2)*(x-x_ref)**1 &
+ 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')
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')
if (dNeq(p1%at(x_ref+x),p1%at(x_ref-x),1e-10_pREAL)) error stop 'polynomials: eval(quadratic)'
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^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 polynomials_selfTest
2022-01-31 19:35:15 +05:30
end module polynomials