Merge branch 'order4-polynomial' into 'development'

allow dependent variables to be of polynomial order up to 4

See merge request damask/DAMASK!572

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), &
+  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,3) * m(1,2) - m(2,3) * m(1,1), &
-               m(1, 2)**2]
+               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), &
+#ifndef __INTEL_LLVM_COMPILER
-               v(2,2) + m(2, 3) * w(1), &
+  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), &
+#ifndef __INTEL_LLVM_COMPILER
-                 v(2,2) + m(2, 3) * w(2), &
+    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])

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 &
+#ifndef __INTEL_LLVM_COMPILER
-                                            + dotState * Delta_t
+  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 &
+#ifndef __INTEL_LLVM_COMPILER
-                                            + dotState * Delta_t
+  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 &
+#ifndef __INTEL_LLVM_COMPILER
-               + A(n,stage) * plastic_RKdotState(1:sizeDotState,n)
+      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 &
+#ifndef __INTEL_LLVM_COMPILER
-                                              + dotState * Delta_t
+    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 &
+#ifndef __INTEL_LLVM_COMPILER
-                                            + dotState * Delta_t
+  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

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)
+  y = self%coef(ubound(self%coef,1))
-  do i = 1, ubound(self%coef,1)
+  do o = ubound(self%coef,1)-1, 0, -1
-    y = y + self%coef(i) * (x-self%x_ref)**i
+#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), dimension(5) :: coef
-  real(pReal) :: x_ref, x
+  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)
+  do i = 1, size(coef_s)
-  write(coef_s(2),*) coef(2)
+    write(coef_s(i),*) coef(i)
-  write(coef_s(3),*) coef(3)
+  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)
+  YAML_s = 'C: 0.0'//IO_EOL//&
-  if (dNeq(p1%at(x_ref+x),-p1%at(x_ref-x),1.0e-10_pReal))         error stop 'polynomials: eval(odd)'
+           'C,T: '//trim(adjustl(coef_s(2)))//IO_EOL//&
-  if (dNeq(p1%der1_at(x),p1%der1_at(5.0_pReal*x),1.0e-10_pReal))  error stop 'polynomials: eval_der(odd)'
+           '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)
+  YAML_s = 'C: 0.0'//IO_EOL//&
-  if (dNeq(p1%at(x_ref+x),p1%at(x_ref-x),1e-10_pReal))            error stop 'polynomials: eval(even)'
+           'C,T^2: '//trim(adjustl(coef_s(3)))//IO_EOL//&
-  if (dNeq(p1%der1_at(x_ref+x),-p1%der1_at(x_ref-x),1e-10_pReal)) error stop 'polynomials: eval_der(even)'
+           '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