From c66e2336c25d704a36a60ea121c837c34210a388 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 14 Feb 2022 08:58:40 +0100 Subject: [PATCH] some versions of ifort have problems with 'do concurrent' --- src/grid/spectral_utilities.f90 | 35 ++++++++++++++++++++++++++++++-- src/math.f90 | 31 ++++++++-------------------- src/phase_mechanical_plastic.f90 | 3 +-- 3 files changed, 42 insertions(+), 27 deletions(-) diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index f89485c7a..2a675c9b6 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -376,20 +376,32 @@ subroutine utilities_updateGamma(C) !$OMP PARALLEL DO PRIVATE(l,m,n,o,temp33_cmplx,xiDyad_cmplx,A,A_inv,err) do k = cells3Offset+1, cells3Offset+cells3; do j = 1, cells(2); do i = 1, cells1Red if (any([i,j,k] /= 1)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 - do concurrent (l = 1:3, m = 1:3) +#ifndef __INTEL_COMPILER + do concurrent(l = 1:3, m = 1:3) xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k-cells3Offset))*xi1st(m,i,j,k-cells3Offset) end do do concurrent(l = 1:3, m = 1:3) temp33_cmplx(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx) end do +#else + forall(l = 1:3, m = 1:3) & + xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k-cells3Offset))*xi1st(m,i,j,k-cells3Offset) + forall(l = 1:3, m = 1:3) & + temp33_cmplx(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx) +#endif A(1:3,1:3) = temp33_cmplx%re; A(4:6,4:6) = temp33_cmplx%re A(1:3,4:6) = temp33_cmplx%im; A(4:6,1:3) = -temp33_cmplx%im if (abs(math_det33(A(1:3,1:3))) > 1e-16) then call math_invert(A_inv, err, A) temp33_cmplx = cmplx(A_inv(1:3,1:3),A_inv(1:3,4:6),pReal) +#ifndef __INTEL_COMPILER do concurrent(l=1:3, m=1:3, n=1:3, o=1:3) gamma_hat(l,m,n,o,i,j,k-cells3Offset) = temp33_cmplx(l,n) * xiDyad_cmplx(o,m) end do +#else + forall(l=1:3, m=1:3, n=1:3, o=1:3) & + gamma_hat(l,m,n,o,i,j,k-cells3Offset) = temp33_cmplx(l,n) * xiDyad_cmplx(o,m) +#endif end if end if end do; end do; end do @@ -497,23 +509,37 @@ subroutine utilities_fourierGammaConvolution(fieldAim) !$OMP PARALLEL DO PRIVATE(l,m,n,o,temp33_cmplx,xiDyad_cmplx,A,A_inv,err,gamma_hat) do k = 1, cells3; do j = 1, cells(2); do i = 1, cells1Red if (any([i,j,k+cells3Offset] /= 1)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 +#ifndef __INTEL_COMPILER do concurrent(l = 1:3, m = 1:3) xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k))*xi1st(m,i,j,k) end do do concurrent(l = 1:3, m = 1:3) temp33_cmplx(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx) end do +#else + forall(l = 1:3, m = 1:3) & + xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k))*xi1st(m,i,j,k) + forall(l = 1:3, m = 1:3) & + temp33_cmplx(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx) +#endif A(1:3,1:3) = temp33_cmplx%re; A(4:6,4:6) = temp33_cmplx%re A(1:3,4:6) = temp33_cmplx%im; A(4:6,1:3) = -temp33_cmplx%im if (abs(math_det33(A(1:3,1:3))) > 1e-16) then call math_invert(A_inv, err, A) temp33_cmplx = cmplx(A_inv(1:3,1:3),A_inv(1:3,4:6),pReal) +#ifndef __INTEL_COMPILER do concurrent(l=1:3, m=1:3, n=1:3, o=1:3) gamma_hat(l,m,n,o,1,1,1) = temp33_cmplx(l,n)*xiDyad_cmplx(o,m) end do do concurrent(l = 1:3, m = 1:3) temp33_cmplx(l,m) = sum(gamma_hat(l,m,1:3,1:3,1,1,1)*tensorField_fourier(1:3,1:3,i,j,k)) end do +#else + forall(l=1:3, m=1:3, n=1:3, o=1:3) & + gamma_hat(l,m,n,o,1,1,1) = temp33_cmplx(l,n)*xiDyad_cmplx(o,m) + forall(l = 1:3, m = 1:3) & + temp33_cmplx(l,m) = sum(gamma_hat(l,m,1:3,1:3,1,1,1)*tensorField_fourier(1:3,1:3,i,j,k)) +#endif tensorField_fourier(1:3,1:3,i,j,k) = temp33_cmplx else tensorField_fourier(1:3,1:3,i,j,k) = cmplx(0.0_pReal,0.0_pReal,pReal) @@ -524,9 +550,14 @@ subroutine utilities_fourierGammaConvolution(fieldAim) else memoryEfficient !$OMP PARALLEL DO PRIVATE(l,m,temp33_cmplx) do k = 1, cells3; do j = 1, cells(2); do i = 1,cells1Red +#ifndef __INTEL_COMPILER do concurrent(l = 1:3, m = 1:3) temp33_cmplx(l,m) = sum(gamma_hat(l,m,1:3,1:3,i,j,k)*tensorField_fourier(1:3,1:3,i,j,k)) end do +#else + forall(l = 1:3, m = 1:3) & + temp33_cmplx(l,m) = sum(gamma_hat(l,m,1:3,1:3,i,j,k)*tensorField_fourier(1:3,1:3,i,j,k)) +#endif tensorField_fourier(1:3,1:3,i,j,k) = temp33_cmplx end do; end do; end do !$OMP END PARALLEL DO @@ -758,7 +789,7 @@ subroutine utilities_fourierVectorDivergence() scalarField_fourier(1:cells1Red,1:cells(2),1:cells3) = sum(vectorField_fourier(1:3,1:cells1Red,1:cells(2),1:cells3) & - *conjg(-xi1st)) + *conjg(-xi1st),1) end subroutine utilities_fourierVectorDivergence diff --git a/src/math.f90 b/src/math.f90 index dd4690672..75b9ddea2 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -262,9 +262,8 @@ pure function math_identity4th() math_identity4th(i,j,k,l) = 0.5_pReal*(math_I3(i,k)*math_I3(j,l)+math_I3(i,l)*math_I3(j,k)) enddo #else - do i=1,3; do j=1,3; do k=1,3; do l=1,3 + forall(i=1:3, j=1:3, k=1:3, l=1:3) & math_identity4th(i,j,k,l) = 0.5_pReal*(math_I3(i,k)*math_I3(j,l)+math_I3(i,l)*math_I3(j,k)) - enddo; enddo; enddo; enddo #endif end function math_identity4th @@ -338,9 +337,7 @@ pure function math_outer(A,B) math_outer(i,j) = A(i)*B(j) enddo #else - do i=1,size(A,1); do j=1,size(B,1) - math_outer(i,j) = A(i)*B(j) - enddo; enddo + forall(i=1:size(A,1), j=1:size(B,1)) math_outer(i,j) = A(i)*B(j) #endif end function math_outer @@ -387,9 +384,7 @@ pure function math_mul3333xx33(A,B) math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3)) enddo #else - do i=1,3; do j=1,3 - math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3)) - enddo; enddo + forall (i=1:3, j=1:3) math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3)) #endif end function math_mul3333xx33 @@ -411,9 +406,7 @@ pure function math_mul3333xx3333(A,B) math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l)) enddo #else - do i=1,3; do j=1,3; do k=1,3; do l=1,3 - math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l)) - enddo; enddo; enddo; enddo + forall(i=1:3, j=1:3, k=1:3, l=1:3) math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l)) #endif end function math_mul3333xx3333 @@ -752,9 +745,7 @@ pure function math_3333to99(m3333) math_3333to99(i,j) = m3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) enddo #else - do i=1,9; do j=1,9 - math_3333to99(i,j) = m3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) - enddo; enddo + forall(i=1:9, j=1:9) math_3333to99(i,j) = m3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) #endif end function math_3333to99 @@ -775,9 +766,7 @@ pure function math_99to3333(m99) math_99to3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) = m99(i,j) enddo #else - do i=1,9; do j=1,9 - math_99to3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) = m99(i,j) - enddo; enddo + forall(i=1:9, j=1:9) math_99to3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) = m99(i,j) #endif end function math_99to3333 @@ -810,9 +799,7 @@ pure function math_sym3333to66(m3333,weighted) math_sym3333to66(i,j) = w(i)*w(j)*m3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j)) enddo #else - do i=1,6; do j=1,6 - math_sym3333to66(i,j) = w(i)*w(j)*m3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j)) - enddo; enddo + forall(i=1:6, j=1:6) math_sym3333to66(i,j) = w(i)*w(j)*m3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j)) #endif end function math_sym3333to66 @@ -950,9 +937,7 @@ pure function math_3333toVoigt66_stiffness(C) result(C_tilde) C_tilde(i,j) = C(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) end do #else - do i=1,6; do j=1,6 - C_tilde(i,j) = C(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) - end do; end do + forall(i=1:6, j=1:6) C_tilde(i,j) = C(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) #endif end function math_3333toVoigt66_stiffness diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index 72b67ef64..3915c3b2d 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -379,10 +379,9 @@ module function plastic_deltaState(ph, en) result(broken) en logical :: broken - real(pReal), dimension(3,3) :: & + real(pReal), dimension(3,3) :: & Mp integer :: & - myOffset, & mySize