diff --git a/src/math.f90 b/src/math.f90 index 94ad2c48d..25c90ccf4 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -82,7 +82,7 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief initialization of random seed generator and internal checks !-------------------------------------------------------------------------------------------------- -subroutine math_init +subroutine math_init() real(pReal), dimension(4) :: randTest integer :: randSize @@ -1045,24 +1045,34 @@ pure subroutine math_eigh33(w,v,m) w = math_eigvalsh33(m) - v(1:3,2) = [ m(1, 2) * m(2, 3) - m(1, 3) * m(2, 2), & - m(1, 3) * m(1, 2) - m(2, 3) * m(1, 1), & - m(1, 2)**2] + v(1:3,2) = [ m(1,2) * m(2,3) - m(1,3) * m(2,2), & + m(1,3) * m(1,2) - m(2,3) * m(1,1), & + m(1,2)**2] T = maxval(abs(w)) U = max(T, T**2) threshold = sqrt(5.68e-14_pReal * U**2) - v(1:3,1) = [ v(1,2) + m(1, 3) * w(1), & - v(2,2) + m(2, 3) * w(1), & +#ifndef __INTEL_LLVM_COMPILER + v(1:3,1) = [m(1,3)*w(1) + v(1,2), & + m(2,3)*w(1) + v(2,2), & +#else + v(1:3,1) = [IEEE_FMA(m(1,3),w(1),v(1,2)), & + IEEE_FMA(m(2,3),w(1),v(2,2)), & +#endif (m(1,1) - w(1)) * (m(2,2) - w(1)) - v(3,2)] norm = norm2(v(1:3, 1)) fallback1: if (norm < threshold) then call math_eigh(w,v,error,m) else fallback1 v(1:3,1) = v(1:3, 1) / norm - v(1:3,2) = [ v(1,2) + m(1, 3) * w(2), & - v(2,2) + m(2, 3) * w(2), & +#ifndef __INTEL_LLVM_COMPILER + v(1:3,2) = [m(1,3)*w(2) + v(1,2), & + m(2,3)*w(2) + v(2,2), & +#else + v(1:3,2) = [IEEE_FMA(m(1,3),w(2),v(1,2)), & + IEEE_FMA(m(2,3),w(2),v(2,2)), & +#endif (m(1,1) - w(2)) * (m(2,2) - w(2)) - v(3,2)] norm = norm2(v(1:3, 2)) fallback2: if (norm < threshold) then @@ -1300,7 +1310,7 @@ end function math_clip !-------------------------------------------------------------------------------------------------- !> @brief Check correctness of some math functions. !-------------------------------------------------------------------------------------------------- -subroutine selfTest +subroutine selfTest() integer, dimension(2,4) :: & sort_in_ = reshape([+1,+5, +5,+6, -1,-1, +3,-2],[2,4]) diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index b13675e38..fc9e40e02 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -680,8 +680,11 @@ function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result if (any(IEEE_is_NaN(dotState))) return sizeDotState = plasticState(ph)%sizeDotState - plasticState(ph)%state(1:sizeDotState,en) = subState0 & - + dotState * Delta_t +#ifndef __INTEL_LLVM_COMPILER + plasticState(ph)%state(1:sizeDotState,en) = subState0 + dotState*Delta_t +#else + plasticState(ph)%state(1:sizeDotState,en) = IEEE_FMA(dotState,Delta_t,subState0) +#endif broken = plastic_deltaState(ph,en) if(broken) return @@ -720,8 +723,11 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en sizeDotState = plasticState(ph)%sizeDotState r = - dotState * 0.5_pReal * Delta_t - plasticState(ph)%state(1:sizeDotState,en) = subState0 & - + dotState * Delta_t +#ifndef __INTEL_LLVM_COMPILER + plasticState(ph)%state(1:sizeDotState,en) = subState0 + dotState*Delta_t +#else + plasticState(ph)%state(1:sizeDotState,en) = IEEE_FMA(dotState,Delta_t,subState0) +#endif broken = plastic_deltaState(ph,en) if(broken) return @@ -842,12 +848,18 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en,A,B,C,DB) dotState = A(1,stage) * plastic_RKdotState(1:sizeDotState,1) do n = 2, stage - dotState = dotState & - + A(n,stage) * plastic_RKdotState(1:sizeDotState,n) +#ifndef __INTEL_LLVM_COMPILER + dotState = dotState + A(n,stage)*plastic_RKdotState(1:sizeDotState,n) +#else + dotState = IEEE_FMA(A(n,stage),plastic_RKdotState(1:sizeDotState,n),dotState) +#endif enddo - plasticState(ph)%state(1:sizeDotState,en) = subState0 & - + dotState * Delta_t +#ifndef __INTEL_LLVM_COMPILER + plasticState(ph)%state(1:sizeDotState,en) = subState0 + dotState*Delta_t +#else + plasticState(ph)%state(1:sizeDotState,en) = IEEE_FMA(dotState,Delta_t,subState0) +#endif broken = integrateStress(F_0+(F-F_0)*Delta_t*C(stage),subFp0,subFi0,Delta_t*C(stage), ph,en) if(broken) exit @@ -861,8 +873,11 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en,A,B,C,DB) plastic_RKdotState(1:sizeDotState,size(B)) = dotState dotState = matmul(plastic_RKdotState,B) - plasticState(ph)%state(1:sizeDotState,en) = subState0 & - + dotState * Delta_t +#ifndef __INTEL_LLVM_COMPILER + plasticState(ph)%state(1:sizeDotState,en) = subState0 + dotState*Delta_t +#else + plasticState(ph)%state(1:sizeDotState,en) = IEEE_FMA(dotState,Delta_t,subState0) +#endif if(present(DB)) & broken = .not. converged(matmul(plastic_RKdotState(1:sizeDotState,1:size(DB)),DB) * Delta_t, & @@ -1146,12 +1161,18 @@ module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF) else lhs_3333 = 0.0_pReal; rhs_3333 = 0.0_pReal do o=1,3; do p=1,3 +#ifndef __INTEL_LLVM_COMPILER lhs_3333(1:3,1:3,o,p) = lhs_3333(1:3,1:3,o,p) & + matmul(invSubFi0,dLidFi(1:3,1:3,o,p)) * Delta_t lhs_3333(1:3,o,1:3,p) = lhs_3333(1:3,o,1:3,p) & + invFi*invFi(p,o) rhs_3333(1:3,1:3,o,p) = rhs_3333(1:3,1:3,o,p) & - matmul(invSubFi0,dLidS(1:3,1:3,o,p)) * Delta_t +#else + lhs_3333(1:3,1:3,o,p) = IEEE_FMA(matmul(invSubFi0,dLidFi(1:3,1:3,o,p)),Delta_t,lhs_3333(1:3,1:3,o,p)) + lhs_3333(1:3,o,1:3,p) = IEEE_FMA(invFi,invFi(p,o),lhs_3333(1:3,o,1:3,p)) + rhs_3333(1:3,1:3,o,p) = IEEE_FMA(matmul(invSubFi0,dLidS(1:3,1:3,o,p)),-Delta_t,rhs_3333(1:3,1:3,o,p)) +#endif enddo; enddo call math_invert(temp_99,error,math_3333to99(lhs_3333)) if (error) then @@ -1180,8 +1201,12 @@ module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF) temp_3333(1:3,1:3,p,o) = matmul(matmul(temp_33_2,dLpdS(1:3,1:3,p,o)), invFi) & + matmul(temp_33_3,dLidS(1:3,1:3,p,o)) enddo; enddo +#ifndef __INTEL_LLVM_COMPILER lhs_3333 = math_mul3333xx3333(dSdFe,temp_3333) * Delta_t & + math_mul3333xx3333(dSdFi,dFidS) +#else + lhs_3333 = IEEE_FMA(math_mul3333xx3333(dSdFe,temp_3333),Delta_t,math_mul3333xx3333(dSdFi,dFidS)) +#endif call math_invert(temp_99,error,math_eye(9)+math_3333to99(lhs_3333)) if (error) then diff --git a/src/polynomials.f90 b/src/polynomials.f90 index 155f54bf8..c883528ca 100644 --- a/src/polynomials.f90 +++ b/src/polynomials.f90 @@ -13,10 +13,9 @@ module polynomials type, public :: tPolynomial real(pReal), dimension(:), allocatable :: coef - real(pReal) :: x_ref + real(pReal) :: x_ref = huge(0.0_pReal) contains procedure, public :: at => eval - procedure, public :: der1_at => eval_der1 end type tPolynomial interface polynomial @@ -46,14 +45,14 @@ end subroutine polynomials_init !-------------------------------------------------------------------------------------------------- !> @brief Initialize a Polynomial from Coefficients. !-------------------------------------------------------------------------------------------------- -function polynomial_from_coef(coef,x_ref) result(p) +pure function polynomial_from_coef(coef,x_ref) result(p) - real(pReal), dimension(:), intent(in) :: coef + real(pReal), dimension(0:), intent(in) :: coef real(pReal), intent(in) :: x_ref type(tPolynomial) :: p - allocate(p%coef(0:size(coef)-1),source=coef) ! should be zero based + p%coef = coef p%x_ref = x_ref end function polynomial_from_coef @@ -70,6 +69,8 @@ function polynomial_from_dict(dict,y,x) result(p) real(pReal), dimension(:), allocatable :: coef real(pReal) :: x_ref + integer :: i, o + character(len=1) :: o_s allocate(coef(1),source=dict%get_asFloat(y)) @@ -77,12 +78,14 @@ function polynomial_from_dict(dict,y,x) result(p) if (dict%contains(y//','//x)) then x_ref = dict%get_asFloat(x//'_ref') coef = [coef,dict%get_asFloat(y//','//x)] - if (dict%contains(y//','//x//'^2')) then - coef = [coef,dict%get_asFloat(y//','//x//'^2')] - end if - else - x_ref = huge(0.0_pReal) ! Simplify debugging end if + do o = 2,4 + write(o_s,'(I0.0)') o + if (dict%contains(y//','//x//'^'//o_s)) then + x_ref = dict%get_asFloat(x//'_ref') + coef = [coef,[(0.0_pReal,i=size(coef),o-1)],dict%get_asFloat(y//','//x//'^'//o_s)] + end if + end do p = Polynomial(coef,x_ref) @@ -91,6 +94,7 @@ end function polynomial_from_dict !-------------------------------------------------------------------------------------------------- !> @brief Evaluate a Polynomial. +!> @details https://nvlpubs.nist.gov/nistpubs/jres/71b/jresv71bn1p11_a1b.pdf (eq. 1.2) !-------------------------------------------------------------------------------------------------- pure function eval(self,x) result(y) @@ -98,49 +102,35 @@ pure function eval(self,x) result(y) real(pReal), intent(in) :: x real(pReal) :: y - integer :: i + integer :: o - y = self%coef(0) - do i = 1, ubound(self%coef,1) - y = y + self%coef(i) * (x-self%x_ref)**i + y = self%coef(ubound(self%coef,1)) + do o = ubound(self%coef,1)-1, 0, -1 +#ifndef __INTEL_LLVM_COMPILER + y = y*(x-self%x_ref) +self%coef(o) +#else + y = IEEE_FMA(y,x-self%x_ref,self%coef(o)) +#endif enddo end function eval -!-------------------------------------------------------------------------------------------------- -!> @brief Evaluate a first derivative of Polynomial. -!-------------------------------------------------------------------------------------------------- -pure function eval_der1(self,x) result(y) - - class(tPolynomial), intent(in) :: self - real(pReal), intent(in) :: x - real(pReal) :: y - - integer :: i - - - y = 0.0_pReal - do i = 1, ubound(self%coef,1) - y = y + real(i,pReal)*self%coef(i) * (x-self%x_ref)**(i-1) - enddo - -end function eval_der1 - - !-------------------------------------------------------------------------------------------------- !> @brief Check correctness of polynomical functionality. !-------------------------------------------------------------------------------------------------- -subroutine selfTest +subroutine selfTest() type(tPolynomial) :: p1, p2 - real(pReal), dimension(3) :: coef - real(pReal) :: x_ref, x + real(pReal), dimension(5) :: coef + integer :: i + real(pReal) :: x_ref, x, y class(tNode), pointer :: dict - character(len=pStringLen), dimension(3) :: coef_s + character(len=pStringLen), dimension(size(coef)) :: coef_s character(len=pStringLen) :: x_ref_s, x_s, YAML_s + call random_number(coef) call random_number(x_ref) call random_number(x) @@ -149,29 +139,56 @@ subroutine selfTest x_ref = x_ref*10_pReal -0.5_pReal x = x*10_pReal -0.5_pReal + p1 = polynomial([coef(1)],x_ref) + if (dNeq(p1%at(x),coef(1))) error stop 'polynomial: eval(constant)' + p1 = polynomial(coef,x_ref) if (dNeq(p1%at(x_ref),coef(1))) error stop 'polynomial: @ref' - write(coef_s(1),*) coef(1) - write(coef_s(2),*) coef(2) - write(coef_s(3),*) coef(3) + do i = 1, size(coef_s) + write(coef_s(i),*) coef(i) + end do 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//& + 'C,T^3: '//trim(adjustl(coef_s(4)))//IO_EOL//& + 'C,T^4: '//trim(adjustl(coef_s(5)))//IO_EOL//& 'T_ref: '//trim(adjustl(x_ref_s))//IO_EOL Dict => YAML_parse_str(trim(YAML_s)) p2 = polynomial(dict%asDict(),'C','T') - if (dNeq(p1%at(x),p2%at(x),1.0e-10_pReal)) error stop 'polynomials: init' + if (dNeq(p1%at(x),p2%at(x),1.0e-6_pReal)) error stop 'polynomials: init' + 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)' - p1 = polynomial(coef*[0.0_pReal,1.0_pReal,0.0_pReal],x_ref) - if (dNeq(p1%at(x_ref+x),-p1%at(x_ref-x),1.0e-10_pReal)) error stop 'polynomials: eval(odd)' - if (dNeq(p1%der1_at(x),p1%der1_at(5.0_pReal*x),1.0e-10_pReal)) error stop 'polynomials: eval_der(odd)' + 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(trim(YAML_s)) + p1 = polynomial(dict%asDict(),'C','T') + if (dNeq(p1%at(x_ref+x),-p1%at(x_ref-x),1.0e-10_pReal)) error stop 'polynomials: eval(linear)' - p1 = polynomial(coef*[0.0_pReal,0.0_pReal,1.0_pReal],x_ref) - if (dNeq(p1%at(x_ref+x),p1%at(x_ref-x),1e-10_pReal)) error stop 'polynomials: eval(even)' - if (dNeq(p1%der1_at(x_ref+x),-p1%der1_at(x_ref-x),1e-10_pReal)) error stop 'polynomials: eval_der(even)' + 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(trim(YAML_s)) + p1 = polynomial(dict%asDict(),'C','T') + if (dNeq(p1%at(x_ref+x),p1%at(x_ref-x),1e-10_pReal)) error stop 'polynomials: eval(quadratic)' + + YAML_s = 'Y: '//trim(adjustl(coef_s(1)))//IO_EOL//& + 'Y,X^3: '//trim(adjustl(coef_s(2)))//IO_EOL//& + 'X_ref: '//trim(adjustl(x_ref_s))//IO_EOL + Dict => YAML_parse_str(trim(YAML_s)) + p1 = polynomial(dict%asDict(),'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)' + + YAML_s = 'Y: '//trim(adjustl(coef_s(1)))//IO_EOL//& + 'Y,X^4: '//trim(adjustl(coef_s(2)))//IO_EOL//& + 'X_ref: '//trim(adjustl(x_ref_s))//IO_EOL + Dict => YAML_parse_str(trim(YAML_s)) + p1 = polynomial(dict%asDict(),'Y','X') + if (dNeq(p1%at(x_ref+x),p1%at(x_ref-x),1.0e-6_pReal)) error stop 'polynomials: eval(quartic)' end subroutine selfTest