From ae3de821b41cdb3b9a0588b5c305bbc7b5b1b835 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 7 Mar 2019 07:01:09 +0100 Subject: [PATCH] bugfix: synchronizatopm of dPdF_min and dPdF_max was not correct -before: using componenwise min/max among different processors -now: identify the processor that holds the minimum/maximum of the norm --- src/spectral_utilities.f90 | 545 ++++++++++++++++++------------------- 1 file changed, 269 insertions(+), 276 deletions(-) diff --git a/src/spectral_utilities.f90 b/src/spectral_utilities.f90 index 877c9aed7..dd753dc7c 100644 --- a/src/spectral_utilities.f90 +++ b/src/spectral_utilities.f90 @@ -584,28 +584,27 @@ end subroutine utilities_fourierGammaConvolution !> @brief doing convolution DamageGreenOp_hat * field_real !-------------------------------------------------------------------------------------------------- subroutine utilities_fourierGreenConvolution(D_ref, mobility_ref, deltaT) - - use math, only: & - math_mul33x3, & - PI - use mesh, only: & - grid, & - grid3 - - implicit none - real(pReal), dimension(3,3), intent(in) :: D_ref !< desired average value of the field after convolution - real(pReal), intent(in) :: mobility_ref, deltaT !< desired average value of the field after convolution - complex(pReal) :: GreenOp_hat - integer(pInt) :: i, j, k - + use math, only: & + math_mul33x3, & + PI + use mesh, only: & + grid, & + grid3 + + implicit none + real(pReal), dimension(3,3), intent(in) :: D_ref + real(pReal), intent(in) :: mobility_ref, deltaT + complex(pReal) :: GreenOp_hat + integer(pInt) :: i, j, k + !-------------------------------------------------------------------------------------------------- ! do the actual spectral method calculation - do k = 1_pInt, grid3; do j = 1_pInt, grid(2) ;do i = 1_pInt, grid1Red - GreenOp_hat = cmplx(1.0_pReal,0.0_pReal,pReal)/ & - (cmplx(mobility_ref,0.0_pReal,pReal) + cmplx(deltaT,0.0_pReal)*& - sum(conjg(xi1st(1:3,i,j,k))* matmul(cmplx(D_ref,0.0_pReal),xi1st(1:3,i,j,k)))) ! why not use dot_product - scalarField_fourier(i,j,k) = scalarField_fourier(i,j,k)*GreenOp_hat - enddo; enddo; enddo + do k = 1_pInt, grid3; do j = 1_pInt, grid(2) ;do i = 1_pInt, grid1Red + GreenOp_hat = cmplx(1.0_pReal,0.0_pReal,pReal)/ & + (cmplx(mobility_ref,0.0_pReal,pReal) + cmplx(deltaT,0.0_pReal)*& + sum(conjg(xi1st(1:3,i,j,k))* matmul(cmplx(D_ref,0.0_pReal),xi1st(1:3,i,j,k)))) ! why not use dot_product + scalarField_fourier(i,j,k) = scalarField_fourier(i,j,k)*GreenOp_hat + enddo; enddo; enddo end subroutine utilities_fourierGreenConvolution @@ -614,47 +613,47 @@ end subroutine utilities_fourierGreenConvolution !> @brief calculate root mean square of divergence of field_fourier !-------------------------------------------------------------------------------------------------- real(pReal) function utilities_divergenceRMS() - use IO, only: & - IO_error - use mesh, only: & - geomSize, & - grid, & - grid3 + use IO, only: & + IO_error + use mesh, only: & + geomSize, & + grid, & + grid3 - implicit none - integer(pInt) :: i, j, k, ierr - complex(pReal), dimension(3) :: rescaledGeom + implicit none + integer(pInt) :: i, j, k, ierr + complex(pReal), dimension(3) :: rescaledGeom - write(6,'(/,a)') ' ... calculating divergence ................................................' - flush(6) + write(6,'(/,a)') ' ... calculating divergence ................................................' + flush(6) - rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal) + rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal) !-------------------------------------------------------------------------------------------------- ! calculating RMS divergence criterion in Fourier space - utilities_divergenceRMS = 0.0_pReal - do k = 1_pInt, grid3; do j = 1_pInt, grid(2) - do i = 2_pInt, grid1Red -1_pInt ! Has somewhere a conj. complex counterpart. Therefore count it twice. - utilities_divergenceRMS = utilities_divergenceRMS & - + 2.0_pReal*(sum (real(matmul(tensorField_fourier(1:3,1:3,i,j,k),& ! (sqrt(real(a)**2 + aimag(a)**2))**2 = real(a)**2 + aimag(a)**2. do not take square root and square again - conjg(-xi1st(1:3,i,j,k))*rescaledGeom))**2.0_pReal)& ! --> sum squared L_2 norm of vector - +sum(aimag(matmul(tensorField_fourier(1:3,1:3,i,j,k),& - conjg(-xi1st(1:3,i,j,k))*rescaledGeom))**2.0_pReal)) - enddo - utilities_divergenceRMS = utilities_divergenceRMS & ! these two layers (DC and Nyquist) do not have a conjugate complex counterpart (if grid(1) /= 1) - + sum( real(matmul(tensorField_fourier(1:3,1:3,1 ,j,k), & - conjg(-xi1st(1:3,1,j,k))*rescaledGeom))**2.0_pReal) & - + sum(aimag(matmul(tensorField_fourier(1:3,1:3,1 ,j,k), & - conjg(-xi1st(1:3,1,j,k))*rescaledGeom))**2.0_pReal) & - + sum( real(matmul(tensorField_fourier(1:3,1:3,grid1Red,j,k), & - conjg(-xi1st(1:3,grid1Red,j,k))*rescaledGeom))**2.0_pReal) & - + sum(aimag(matmul(tensorField_fourier(1:3,1:3,grid1Red,j,k), & - conjg(-xi1st(1:3,grid1Red,j,k))*rescaledGeom))**2.0_pReal) - enddo; enddo - if(grid(1) == 1_pInt) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pReal ! counted twice in case of grid(1) == 1 - call MPI_Allreduce(MPI_IN_PLACE,utilities_divergenceRMS,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='utilities_divergenceRMS') - utilities_divergenceRMS = sqrt(utilities_divergenceRMS) * wgt ! RMS in real space calculated with Parsevals theorem from Fourier space + utilities_divergenceRMS = 0.0_pReal + do k = 1_pInt, grid3; do j = 1_pInt, grid(2) + do i = 2_pInt, grid1Red -1_pInt ! Has somewhere a conj. complex counterpart. Therefore count it twice. + utilities_divergenceRMS = utilities_divergenceRMS & + + 2.0_pReal*(sum (real(matmul(tensorField_fourier(1:3,1:3,i,j,k),& ! (sqrt(real(a)**2 + aimag(a)**2))**2 = real(a)**2 + aimag(a)**2. do not take square root and square again + conjg(-xi1st(1:3,i,j,k))*rescaledGeom))**2.0_pReal)& ! --> sum squared L_2 norm of vector + +sum(aimag(matmul(tensorField_fourier(1:3,1:3,i,j,k),& + conjg(-xi1st(1:3,i,j,k))*rescaledGeom))**2.0_pReal)) + enddo + utilities_divergenceRMS = utilities_divergenceRMS & ! these two layers (DC and Nyquist) do not have a conjugate complex counterpart (if grid(1) /= 1) + + sum( real(matmul(tensorField_fourier(1:3,1:3,1 ,j,k), & + conjg(-xi1st(1:3,1,j,k))*rescaledGeom))**2.0_pReal) & + + sum(aimag(matmul(tensorField_fourier(1:3,1:3,1 ,j,k), & + conjg(-xi1st(1:3,1,j,k))*rescaledGeom))**2.0_pReal) & + + sum( real(matmul(tensorField_fourier(1:3,1:3,grid1Red,j,k), & + conjg(-xi1st(1:3,grid1Red,j,k))*rescaledGeom))**2.0_pReal) & + + sum(aimag(matmul(tensorField_fourier(1:3,1:3,grid1Red,j,k), & + conjg(-xi1st(1:3,grid1Red,j,k))*rescaledGeom))**2.0_pReal) + enddo; enddo + if(grid(1) == 1_pInt) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pReal ! counted twice in case of grid(1) == 1 + call MPI_Allreduce(MPI_IN_PLACE,utilities_divergenceRMS,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='utilities_divergenceRMS') + utilities_divergenceRMS = sqrt(utilities_divergenceRMS) * wgt ! RMS in real space calculated with Parsevals theorem from Fourier space end function utilities_divergenceRMS @@ -664,66 +663,66 @@ end function utilities_divergenceRMS !> @brief calculate max of curl of field_fourier !-------------------------------------------------------------------------------------------------- real(pReal) function utilities_curlRMS() - use IO, only: & - IO_error - use mesh, only: & - geomSize, & - grid, & - grid3 - - implicit none - integer(pInt) :: i, j, k, l, ierr - complex(pReal), dimension(3,3) :: curl_fourier - complex(pReal), dimension(3) :: rescaledGeom - - write(6,'(/,a)') ' ... calculating curl ......................................................' - flush(6) - - rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal) - + use IO, only: & + IO_error + use mesh, only: & + geomSize, & + grid, & + grid3 + + implicit none + integer(pInt) :: i, j, k, l, ierr + complex(pReal), dimension(3,3) :: curl_fourier + complex(pReal), dimension(3) :: rescaledGeom + + write(6,'(/,a)') ' ... calculating curl ......................................................' + flush(6) + + rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal) + !-------------------------------------------------------------------------------------------------- ! calculating max curl criterion in Fourier space - utilities_curlRMS = 0.0_pReal - - do k = 1_pInt, grid3; do j = 1_pInt, grid(2); - do i = 2_pInt, grid1Red - 1_pInt - do l = 1_pInt, 3_pInt - curl_fourier(l,1) = (+tensorField_fourier(l,3,i,j,k)*xi1st(2,i,j,k)*rescaledGeom(2) & - -tensorField_fourier(l,2,i,j,k)*xi1st(3,i,j,k)*rescaledGeom(3)) - curl_fourier(l,2) = (+tensorField_fourier(l,1,i,j,k)*xi1st(3,i,j,k)*rescaledGeom(3) & - -tensorField_fourier(l,3,i,j,k)*xi1st(1,i,j,k)*rescaledGeom(1)) - curl_fourier(l,3) = (+tensorField_fourier(l,2,i,j,k)*xi1st(1,i,j,k)*rescaledGeom(1) & - -tensorField_fourier(l,1,i,j,k)*xi1st(2,i,j,k)*rescaledGeom(2)) - enddo - utilities_curlRMS = utilities_curlRMS & - +2.0_pReal*sum(real(curl_fourier)**2.0_pReal+aimag(curl_fourier)**2.0_pReal) ! Has somewhere a conj. complex counterpart. Therefore count it twice. - enddo - do l = 1_pInt, 3_pInt - curl_fourier = (+tensorField_fourier(l,3,1,j,k)*xi1st(2,1,j,k)*rescaledGeom(2) & - -tensorField_fourier(l,2,1,j,k)*xi1st(3,1,j,k)*rescaledGeom(3)) - curl_fourier = (+tensorField_fourier(l,1,1,j,k)*xi1st(3,1,j,k)*rescaledGeom(3) & - -tensorField_fourier(l,3,1,j,k)*xi1st(1,1,j,k)*rescaledGeom(1)) - curl_fourier = (+tensorField_fourier(l,2,1,j,k)*xi1st(1,1,j,k)*rescaledGeom(1) & - -tensorField_fourier(l,1,1,j,k)*xi1st(2,1,j,k)*rescaledGeom(2)) - enddo - utilities_curlRMS = utilities_curlRMS & - + sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal) ! this layer (DC) does not have a conjugate complex counterpart (if grid(1) /= 1) - do l = 1_pInt, 3_pInt - curl_fourier = (+tensorField_fourier(l,3,grid1Red,j,k)*xi1st(2,grid1Red,j,k)*rescaledGeom(2) & - -tensorField_fourier(l,2,grid1Red,j,k)*xi1st(3,grid1Red,j,k)*rescaledGeom(3)) - curl_fourier = (+tensorField_fourier(l,1,grid1Red,j,k)*xi1st(3,grid1Red,j,k)*rescaledGeom(3) & - -tensorField_fourier(l,3,grid1Red,j,k)*xi1st(1,grid1Red,j,k)*rescaledGeom(1)) - curl_fourier = (+tensorField_fourier(l,2,grid1Red,j,k)*xi1st(1,grid1Red,j,k)*rescaledGeom(1) & - -tensorField_fourier(l,1,grid1Red,j,k)*xi1st(2,grid1Red,j,k)*rescaledGeom(2)) - enddo - utilities_curlRMS = utilities_curlRMS & - + sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal) ! this layer (Nyquist) does not have a conjugate complex counterpart (if grid(1) /= 1) - enddo; enddo - - call MPI_Allreduce(MPI_IN_PLACE,utilities_curlRMS,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='utilities_curlRMS') - utilities_curlRMS = sqrt(utilities_curlRMS) * wgt - if(grid(1) == 1_pInt) utilities_curlRMS = utilities_curlRMS * 0.5_pReal ! counted twice in case of grid(1) == 1 + utilities_curlRMS = 0.0_pReal + + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); + do i = 2_pInt, grid1Red - 1_pInt + do l = 1_pInt, 3_pInt + curl_fourier(l,1) = (+tensorField_fourier(l,3,i,j,k)*xi1st(2,i,j,k)*rescaledGeom(2) & + -tensorField_fourier(l,2,i,j,k)*xi1st(3,i,j,k)*rescaledGeom(3)) + curl_fourier(l,2) = (+tensorField_fourier(l,1,i,j,k)*xi1st(3,i,j,k)*rescaledGeom(3) & + -tensorField_fourier(l,3,i,j,k)*xi1st(1,i,j,k)*rescaledGeom(1)) + curl_fourier(l,3) = (+tensorField_fourier(l,2,i,j,k)*xi1st(1,i,j,k)*rescaledGeom(1) & + -tensorField_fourier(l,1,i,j,k)*xi1st(2,i,j,k)*rescaledGeom(2)) + enddo + utilities_curlRMS = utilities_curlRMS & + +2.0_pReal*sum(real(curl_fourier)**2.0_pReal+aimag(curl_fourier)**2.0_pReal)! Has somewhere a conj. complex counterpart. Therefore count it twice. + enddo + do l = 1_pInt, 3_pInt + curl_fourier = (+tensorField_fourier(l,3,1,j,k)*xi1st(2,1,j,k)*rescaledGeom(2) & + -tensorField_fourier(l,2,1,j,k)*xi1st(3,1,j,k)*rescaledGeom(3)) + curl_fourier = (+tensorField_fourier(l,1,1,j,k)*xi1st(3,1,j,k)*rescaledGeom(3) & + -tensorField_fourier(l,3,1,j,k)*xi1st(1,1,j,k)*rescaledGeom(1)) + curl_fourier = (+tensorField_fourier(l,2,1,j,k)*xi1st(1,1,j,k)*rescaledGeom(1) & + -tensorField_fourier(l,1,1,j,k)*xi1st(2,1,j,k)*rescaledGeom(2)) + enddo + utilities_curlRMS = utilities_curlRMS & + + sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal) ! this layer (DC) does not have a conjugate complex counterpart (if grid(1) /= 1) + do l = 1_pInt, 3_pInt + curl_fourier = (+tensorField_fourier(l,3,grid1Red,j,k)*xi1st(2,grid1Red,j,k)*rescaledGeom(2) & + -tensorField_fourier(l,2,grid1Red,j,k)*xi1st(3,grid1Red,j,k)*rescaledGeom(3)) + curl_fourier = (+tensorField_fourier(l,1,grid1Red,j,k)*xi1st(3,grid1Red,j,k)*rescaledGeom(3) & + -tensorField_fourier(l,3,grid1Red,j,k)*xi1st(1,grid1Red,j,k)*rescaledGeom(1)) + curl_fourier = (+tensorField_fourier(l,2,grid1Red,j,k)*xi1st(1,grid1Red,j,k)*rescaledGeom(1) & + -tensorField_fourier(l,1,grid1Red,j,k)*xi1st(2,grid1Red,j,k)*rescaledGeom(2)) + enddo + utilities_curlRMS = utilities_curlRMS & + + sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal) ! this layer (Nyquist) does not have a conjugate complex counterpart (if grid(1) /= 1) + enddo; enddo + + call MPI_Allreduce(MPI_IN_PLACE,utilities_curlRMS,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='utilities_curlRMS') + utilities_curlRMS = sqrt(utilities_curlRMS) * wgt + if(grid(1) == 1_pInt) utilities_curlRMS = utilities_curlRMS * 0.5_pReal ! counted twice in case of grid(1) == 1 end function utilities_curlRMS @@ -817,9 +816,6 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C) write(6,trim(formatString),advance='no') ' S (load) ', transpose(s_reduced) if(errmatinv) call IO_error(error_ID=400_pInt,ext_msg='utilities_maskedCompliance') endif - deallocate(c_reduced) - deallocate(s_reduced) - deallocate(sTimesC) else temp99_real = 0.0_pReal endif @@ -837,17 +833,17 @@ end function utilities_maskedCompliance !> @brief calculate scalar gradient in fourier field !-------------------------------------------------------------------------------------------------- subroutine utilities_fourierScalarGradient() - use mesh, only: & - grid3, & - grid - - implicit none - integer(pInt) :: i, j, k - - vectorField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal) - forall(k = 1_pInt:grid3, j = 1_pInt:grid(2), i = 1_pInt:grid1Red) & - vectorField_fourier(1:3,i,j,k) = scalarField_fourier(i,j,k)*xi1st(1:3,i,j,k) - + use mesh, only: & + grid3, & + grid + + implicit none + integer(pInt) :: i, j, k + + vectorField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal) + forall(k = 1_pInt:grid3, j = 1_pInt:grid(2), i = 1_pInt:grid1Red) & + vectorField_fourier(1:3,i,j,k) = scalarField_fourier(i,j,k)*xi1st(1:3,i,j,k) + end subroutine utilities_fourierScalarGradient @@ -855,18 +851,18 @@ end subroutine utilities_fourierScalarGradient !> @brief calculate vector divergence in fourier field !-------------------------------------------------------------------------------------------------- subroutine utilities_fourierVectorDivergence() - use mesh, only: & - grid3, & - grid - - implicit none - integer(pInt) :: i, j, k - - scalarField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal) - forall(k = 1_pInt:grid3, j = 1_pInt:grid(2), i = 1_pInt:grid1Red) & - scalarField_fourier(i,j,k) = scalarField_fourier(i,j,k) + & - sum(vectorField_fourier(1:3,i,j,k)*conjg(-xi1st(1:3,i,j,k))) - + use mesh, only: & + grid3, & + grid + + implicit none + integer(pInt) :: i, j, k + + scalarField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal) + forall(k = 1_pInt:grid3, j = 1_pInt:grid(2), i = 1_pInt:grid1Red) & + scalarField_fourier(i,j,k) = scalarField_fourier(i,j,k) + & + sum(vectorField_fourier(1:3,i,j,k)*conjg(-xi1st(1:3,i,j,k))) + end subroutine utilities_fourierVectorDivergence @@ -874,19 +870,20 @@ end subroutine utilities_fourierVectorDivergence !> @brief calculate vector gradient in fourier field !-------------------------------------------------------------------------------------------------- subroutine utilities_fourierVectorGradient() - use mesh, only: & - grid3, & - grid + use mesh, only: & + grid3, & + grid + + implicit none + integer(pInt) :: i, j, k, m, n + + tensorField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal) + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid1Red + do m = 1_pInt, 3_pInt; do n = 1_pInt, 3_pInt + tensorField_fourier(m,n,i,j,k) = vectorField_fourier(m,i,j,k)*xi1st(n,i,j,k) + enddo; enddo + enddo; enddo; enddo - implicit none - integer(pInt) :: i, j, k, m, n - - tensorField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal) - do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid1Red - do m = 1_pInt, 3_pInt; do n = 1_pInt, 3_pInt - tensorField_fourier(m,n,i,j,k) = vectorField_fourier(m,i,j,k)*xi1st(n,i,j,k) - enddo; enddo - enddo; enddo; enddo end subroutine utilities_fourierVectorGradient @@ -894,21 +891,22 @@ end subroutine utilities_fourierVectorGradient !> @brief calculate tensor divergence in fourier field !-------------------------------------------------------------------------------------------------- subroutine utilities_fourierTensorDivergence() - use mesh, only: & - grid3, & - grid + use mesh, only: & + grid3, & + grid + + implicit none + integer(pInt) :: i, j, k, m, n + + vectorField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal) + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid1Red + do m = 1_pInt, 3_pInt; do n = 1_pInt, 3_pInt + vectorField_fourier(m,i,j,k) = & + vectorField_fourier(m,i,j,k) + & + tensorField_fourier(m,n,i,j,k)*conjg(-xi1st(n,i,j,k)) + enddo; enddo + enddo; enddo; enddo - implicit none - integer(pInt) :: i, j, k, m, n - - vectorField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal) - do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid1Red - do m = 1_pInt, 3_pInt; do n = 1_pInt, 3_pInt - vectorField_fourier(m,i,j,k) = & - vectorField_fourier(m,i,j,k) + & - tensorField_fourier(m,n,i,j,k)*conjg(-xi1st(n,i,j,k)) - enddo; enddo - enddo; enddo; enddo end subroutine utilities_fourierTensorDivergence @@ -917,99 +915,93 @@ end subroutine utilities_fourierTensorDivergence !-------------------------------------------------------------------------------------------------- subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& F,timeinc,rotation_BC) - use IO, only: & - IO_error - use debug, only: & - debug_reset, & - debug_info - use math, only: & - math_rotate_forward33, & - math_det33 - use mesh, only: & - grid,& - grid3 - use homogenization, only: & - materialpoint_F, & - materialpoint_P, & - materialpoint_dPdF, & - materialpoint_stressAndItsTangent - - implicit none - real(pReal),intent(out), dimension(3,3,3,3) :: C_volAvg, C_minmaxAvg !< average stiffness - real(pReal),intent(out), dimension(3,3) :: P_av !< average PK stress - real(pReal),intent(out), dimension(3,3,grid(1),grid(2),grid3) :: P !< PK stress - - real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: F !< deformation gradient target !< previous deformation gradient - real(pReal), intent(in) :: timeinc !< loading time - real(pReal), intent(in), dimension(3,3) :: rotation_BC !< rotation of load frame - + use IO, only: & + IO_error + use numerics, only: & + worldrank + use debug, only: & + debug_reset, & + debug_info + use math, only: & + math_rotate_forward33, & + math_det33 + use mesh, only: & + grid,& + grid3 + use homogenization, only: & + materialpoint_F, & + materialpoint_P, & + materialpoint_dPdF, & + materialpoint_stressAndItsTangent - integer(pInt) :: & - j,k,ierr - real(pReal), dimension(3,3,3,3) :: max_dPdF, min_dPdF - real(pReal) :: max_dPdF_norm, min_dPdF_norm, defgradDetMin, defgradDetMax, defgradDet - - write(6,'(/,a)') ' ... evaluating constitutive response ......................................' - flush(6) + implicit none + real(pReal),intent(out), dimension(3,3,3,3) :: C_volAvg, C_minmaxAvg !< average stiffness + real(pReal),intent(out), dimension(3,3) :: P_av !< average PK stress + real(pReal),intent(out), dimension(3,3,grid(1),grid(2),grid3) :: P !< PK stress + real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: F !< deformation gradient target + real(pReal), intent(in) :: timeinc !< loading time + real(pReal), intent(in), dimension(3,3) :: rotation_BC !< rotation of load frame - materialpoint_F = reshape(F,[3,3,1,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field + + integer(pInt) :: & + i,ierr + real(pReal), dimension(3,3,3,3) :: dPdF_max, dPdF_min + real(pReal) :: dPdF_norm_max, dPdF_norm_min + real(pReal), dimension(2) :: valueAndRank !< pair of min/max norm of dPdF to synchronize min/max of dPdF -!-------------------------------------------------------------------------------------------------- -! calculate bounds of det(F) and report - if(debugGeneral) then - defgradDetMax = -huge(1.0_pReal) - defgradDetMin = +huge(1.0_pReal) - do j = 1_pInt, product(grid(1:2))*grid3 - defgradDet = math_det33(materialpoint_F(1:3,1:3,1,j)) - defgradDetMax = max(defgradDetMax,defgradDet) - defgradDetMin = min(defgradDetMin,defgradDet) - end do - - write(6,'(a,1x,es11.4)') ' max determinant of deformation =', defgradDetMax - write(6,'(a,1x,es11.4)') ' min determinant of deformation =', defgradDetMin - flush(6) - endif - - call debug_reset() ! this has no effect on rank >0 - call materialpoint_stressAndItsTangent(.true.,timeinc) ! calculate P field - - P = reshape(materialpoint_P, [3,3,grid(1),grid(2),grid3]) - P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt ! average of P - call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - if (debugRotation) & - write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress (lab) / MPa =',& - transpose(P_av)*1.e-6_pReal - P_av = math_rotate_forward33(P_av,rotation_BC) - write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',& - transpose(P_av)*1.e-6_pReal - flush(6) - - max_dPdF = 0.0_pReal - max_dPdF_norm = 0.0_pReal - min_dPdF = huge(1.0_pReal) - min_dPdF_norm = huge(1.0_pReal) - do k = 1_pInt, product(grid(1:2))*grid3 - if (max_dPdF_norm < sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,k)**2.0_pReal)) then - max_dPdF = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,k) - max_dPdF_norm = sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,k)**2.0_pReal) - endif - if (min_dPdF_norm > sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,k)**2.0_pReal)) then - min_dPdF = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,k) - min_dPdF_norm = sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,k)**2.0_pReal) - endif - end do - - call MPI_Allreduce(MPI_IN_PLACE,max_dPdF,81,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr) - if (ierr /= 0_pInt) call IO_error(894_pInt, ext_msg='MPI_Allreduce max') - call MPI_Allreduce(MPI_IN_PLACE,min_dPdF,81,MPI_DOUBLE,MPI_MIN,PETSC_COMM_WORLD,ierr) - if (ierr /= 0_pInt) call IO_error(894_pInt, ext_msg='MPI_Allreduce min') - - C_minmaxAvg = 0.5_pReal*(max_dPdF + min_dPdF) - - C_volAvg = sum(sum(materialpoint_dPdF,dim=6),dim=5) * wgt - call MPI_Allreduce(MPI_IN_PLACE,C_volAvg,81,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - - call debug_info() ! this has no effect on rank >0 + write(6,'(/,a)') ' ... evaluating constitutive response ......................................' + flush(6) + + materialpoint_F = reshape(F,[3,3,1,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field + + call debug_reset() ! this has no effect on rank >0 + call materialpoint_stressAndItsTangent(.true.,timeinc) ! calculate P field + + P = reshape(materialpoint_P, [3,3,grid(1),grid(2),grid3]) + P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt ! average of P + call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + if (debugRotation) & + write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress (lab) / MPa =',& + transpose(P_av)*1.e-6_pReal + P_av = math_rotate_forward33(P_av,rotation_BC) + write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',& + transpose(P_av)*1.e-6_pReal + flush(6) + + dPdF_max = 0.0_pReal + dPdF_norm_max = 0.0_pReal + dPdF_min = huge(1.0_pReal) + dPdF_norm_min = huge(1.0_pReal) + do i = 1_pInt, product(grid(1:2))*grid3 + if (dPdF_norm_max < sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal)) then + dPdF_max = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,i) + dPdF_norm_max = sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal) + endif + if (dPdF_norm_min > sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal)) then + dPdF_min = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,i) + dPdF_norm_min = sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal) + endif + end do + + valueAndRank = [dPdF_norm_max,real(worldrank,pReal)] + call MPI_Allreduce(MPI_IN_PLACE,valueAndRank,1, MPI_2DOUBLE_PRECISION, MPI_MAXLOC, PETSC_COMM_WORLD, ierr) + if (ierr /= 0_pInt) call IO_error(894_pInt, ext_msg='MPI_Allreduce max') + call MPI_Bcast(dPdF_max,81,MPI_DOUBLE,int(valueAndRank(2)),PETSC_COMM_WORLD, ierr) + if (ierr /= 0_pInt) call IO_error(894_pInt, ext_msg='MPI_Bcast max') + + valueAndRank = [dPdF_norm_min,real(worldrank,pReal)] + call MPI_Allreduce(MPI_IN_PLACE,valueAndRank,1, MPI_2DOUBLE_PRECISION, MPI_MINLOC, PETSC_COMM_WORLD, ierr) + if (ierr /= 0_pInt) call IO_error(894_pInt, ext_msg='MPI_Allreduce min') + call MPI_Bcast(dPdF_min,81,MPI_DOUBLE,int(valueAndRank(2)),PETSC_COMM_WORLD, ierr) + if (ierr /= 0_pInt) call IO_error(894_pInt, ext_msg='MPI_Bcast min') + + C_minmaxAvg = 0.5_pReal*(dPdF_max + dPdF_min) + + C_volAvg = sum(sum(materialpoint_dPdF,dim=6),dim=5) + call MPI_Allreduce(MPI_IN_PLACE,C_volAvg,81,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + C_volAvg = C_volAvg * wgt + + call debug_info() ! this has no effect on rank >0 end subroutine utilities_constitutiveResponse @@ -1018,27 +1010,28 @@ end subroutine utilities_constitutiveResponse !> @brief calculates forward rate, either guessing or just add delta/timeinc !-------------------------------------------------------------------------------------------------- pure function utilities_calculateRate(heterogeneous,field0,field,dt,avRate) - use mesh, only: & - grid3, & - grid - - implicit none - real(pReal), intent(in), dimension(3,3) :: avRate !< homogeneous addon - real(pReal), intent(in) :: & - dt !< timeinc between field0 and field - logical, intent(in) :: & - heterogeneous !< calculate field of rates - real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: & - field0, & !< data of previous step - field !< data of current step - real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: & - utilities_calculateRate - - if (heterogeneous) then - utilities_calculateRate = (field-field0) / dt - else - utilities_calculateRate = spread(spread(spread(avRate,3,grid(1)),4,grid(2)),5,grid3) - endif + use mesh, only: & + grid3, & + grid + + implicit none + real(pReal), intent(in), dimension(3,3) :: & + avRate !< homogeneous addon + real(pReal), intent(in) :: & + dt !< timeinc between field0 and field + logical, intent(in) :: & + heterogeneous !< calculate field of rates + real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: & + field0, & !< data of previous step + field !< data of current step + real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: & + utilities_calculateRate + + if (heterogeneous) then + utilities_calculateRate = (field-field0) / dt + else + utilities_calculateRate = spread(spread(spread(avRate,3,grid(1)),4,grid(2)),5,grid3) + endif end function utilities_calculateRate