From ae3de821b41cdb3b9a0588b5c305bbc7b5b1b835 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 7 Mar 2019 07:01:09 +0100 Subject: [PATCH 1/4] 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 From 21bbba1575275556f249e7ce170d9787243ad686 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 7 Mar 2019 07:09:00 +0100 Subject: [PATCH 2/4] cleaning --- src/numerics.f90 | 58 ++------------------------------------ src/spectral_utilities.f90 | 8 ------ 2 files changed, 2 insertions(+), 64 deletions(-) diff --git a/src/numerics.f90 b/src/numerics.f90 index 9727a04a7..fbc5f52dc 100644 --- a/src/numerics.f90 +++ b/src/numerics.f90 @@ -70,22 +70,12 @@ module numerics err_thermal_tolAbs = 1.0e-2_pReal, & !< absolute tolerance for thermal equilibrium err_thermal_tolRel = 1.0e-6_pReal, & !< relative tolerance for thermal equilibrium err_damage_tolAbs = 1.0e-2_pReal, & !< absolute tolerance for damage evolution - err_damage_tolRel = 1.0e-6_pReal, & !< relative tolerance for damage evolution - err_vacancyflux_tolAbs = 1.0e-8_pReal, & !< absolute tolerance for vacancy transport - err_vacancyflux_tolRel = 1.0e-6_pReal, & !< relative tolerance for vacancy transport - err_porosity_tolAbs = 1.0e-2_pReal, & !< absolute tolerance for porosity evolution - err_porosity_tolRel = 1.0e-6_pReal, & !< relative tolerance for porosity evolution - err_hydrogenflux_tolAbs = 1.0e-8_pReal, & !< absolute tolerance for hydrogen transport - err_hydrogenflux_tolRel = 1.0e-6_pReal, & !< relative tolerance for hydrogen transport - vacancyBoundPenalty = 1.0e+4_pReal, & !< penalty to enforce 0 < Cv < 1 - hydrogenBoundPenalty = 1.0e+4_pReal !< penalty to enforce 0 < Ch < 1 + err_damage_tolRel = 1.0e-6_pReal !< relative tolerance for damage evolution integer(pInt), protected, public :: & itmax = 250_pInt, & !< maximum number of iterations itmin = 1_pInt, & !< minimum number of iterations stagItMax = 10_pInt, & !< max number of field level staggered iterations - maxCutBack = 3_pInt, & !< max number of cut backs - vacancyPolyOrder = 10_pInt, & !< order of polynomial approximation of entropic contribution to vacancy chemical potential - hydrogenPolyOrder = 10_pInt !< order of polynomial approximation of entropic contribution to hydrogen chemical potential + maxCutBack = 3_pInt !< max number of cut backs !-------------------------------------------------------------------------------------------------- ! spectral parameters: @@ -153,11 +143,6 @@ contains ! a sanity check !-------------------------------------------------------------------------------------------------- subroutine numerics_init -#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 - use, intrinsic :: iso_fortran_env, only: & - compiler_version, & - compiler_options -#endif use IO, only: & IO_read, & IO_error, & @@ -191,8 +176,6 @@ subroutine numerics_init call MPI_Comm_size(PETSC_COMM_WORLD,worldsize,ierr);CHKERRQ(ierr) #endif write(6,'(/,a)') ' <<<+- numerics init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() -#include "compilation_info.f90" !$ call GET_ENVIRONMENT_VARIABLE(NAME='DAMASK_NUM_THREADS',VALUE=DAMASK_NumThreadsString,STATUS=gotDAMASK_NUM_THREADS) ! get environment variable DAMASK_NUM_THREADS... !$ if(gotDAMASK_NUM_THREADS /= 0) then ! could not get number of threads, set it to 1 @@ -327,22 +310,6 @@ subroutine numerics_init err_damage_tolabs = IO_floatValue(line,chunkPos,2_pInt) case ('err_damage_tolrel') err_damage_tolrel = IO_floatValue(line,chunkPos,2_pInt) - case ('err_vacancyflux_tolabs') - err_vacancyflux_tolabs = IO_floatValue(line,chunkPos,2_pInt) - case ('err_vacancyflux_tolrel') - err_vacancyflux_tolrel = IO_floatValue(line,chunkPos,2_pInt) - case ('err_porosity_tolabs') - err_porosity_tolabs = IO_floatValue(line,chunkPos,2_pInt) - case ('err_porosity_tolrel') - err_porosity_tolrel = IO_floatValue(line,chunkPos,2_pInt) - case ('err_hydrogenflux_tolabs') - err_hydrogenflux_tolabs = IO_floatValue(line,chunkPos,2_pInt) - case ('err_hydrogenflux_tolrel') - err_hydrogenflux_tolrel = IO_floatValue(line,chunkPos,2_pInt) - case ('vacancyboundpenalty') - vacancyBoundPenalty = IO_floatValue(line,chunkPos,2_pInt) - case ('hydrogenboundpenalty') - hydrogenBoundPenalty = IO_floatValue(line,chunkPos,2_pInt) case ('itmax') itmax = IO_intValue(line,chunkPos,2_pInt) case ('itmin') @@ -351,10 +318,6 @@ subroutine numerics_init maxCutBack = IO_intValue(line,chunkPos,2_pInt) case ('maxstaggerediter') stagItMax = IO_intValue(line,chunkPos,2_pInt) - case ('vacancypolyorder') - vacancyPolyOrder = IO_intValue(line,chunkPos,2_pInt) - case ('hydrogenpolyorder') - hydrogenPolyOrder = IO_intValue(line,chunkPos,2_pInt) !-------------------------------------------------------------------------------------------------- ! spectral parameters @@ -509,22 +472,12 @@ subroutine numerics_init write(6,'(a24,1x,i8)') ' itmin: ',itmin write(6,'(a24,1x,i8)') ' maxCutBack: ',maxCutBack write(6,'(a24,1x,i8)') ' maxStaggeredIter: ',stagItMax - write(6,'(a24,1x,i8)') ' vacancyPolyOrder: ',vacancyPolyOrder - write(6,'(a24,1x,i8)') ' hydrogenPolyOrder: ',hydrogenPolyOrder write(6,'(a24,1x,es8.1)') ' err_struct_tolAbs: ',err_struct_tolAbs write(6,'(a24,1x,es8.1)') ' err_struct_tolRel: ',err_struct_tolRel write(6,'(a24,1x,es8.1)') ' err_thermal_tolabs: ',err_thermal_tolabs write(6,'(a24,1x,es8.1)') ' err_thermal_tolrel: ',err_thermal_tolrel write(6,'(a24,1x,es8.1)') ' err_damage_tolabs: ',err_damage_tolabs write(6,'(a24,1x,es8.1)') ' err_damage_tolrel: ',err_damage_tolrel - write(6,'(a24,1x,es8.1)') ' err_vacancyflux_tolabs: ',err_vacancyflux_tolabs - write(6,'(a24,1x,es8.1)') ' err_vacancyflux_tolrel: ',err_vacancyflux_tolrel - write(6,'(a24,1x,es8.1)') ' err_porosity_tolabs: ',err_porosity_tolabs - write(6,'(a24,1x,es8.1)') ' err_porosity_tolrel: ',err_porosity_tolrel - write(6,'(a24,1x,es8.1)') ' err_hydrogenflux_tolabs:',err_hydrogenflux_tolabs - write(6,'(a24,1x,es8.1)') ' err_hydrogenflux_tolrel:',err_hydrogenflux_tolrel - write(6,'(a24,1x,es8.1)') ' vacancyBoundPenalty: ',vacancyBoundPenalty - write(6,'(a24,1x,es8.1)') ' hydrogenBoundPenalty: ',hydrogenBoundPenalty !-------------------------------------------------------------------------------------------------- ! spectral parameters @@ -608,19 +561,12 @@ subroutine numerics_init if (itmin > itmax .or. itmin < 1_pInt) call IO_error(301_pInt,ext_msg='itmin') if (maxCutBack < 0_pInt) call IO_error(301_pInt,ext_msg='maxCutBack') if (stagItMax < 0_pInt) call IO_error(301_pInt,ext_msg='maxStaggeredIter') - if (vacancyPolyOrder < 0_pInt) call IO_error(301_pInt,ext_msg='vacancyPolyOrder') if (err_struct_tolRel <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_struct_tolRel') if (err_struct_tolAbs <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_struct_tolAbs') if (err_thermal_tolabs <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_thermal_tolabs') if (err_thermal_tolrel <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_thermal_tolrel') if (err_damage_tolabs <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_damage_tolabs') if (err_damage_tolrel <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_damage_tolrel') - if (err_vacancyflux_tolabs <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_vacancyflux_tolabs') - if (err_vacancyflux_tolrel <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_vacancyflux_tolrel') - if (err_porosity_tolabs <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_porosity_tolabs') - if (err_porosity_tolrel <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_porosity_tolrel') - if (err_hydrogenflux_tolabs <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_hydrogenflux_tolabs') - if (err_hydrogenflux_tolrel <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_hydrogenflux_tolrel') #ifdef Spectral if (divergence_correction < 0_pInt .or. & divergence_correction > 2_pInt) call IO_error(301_pInt,ext_msg='divergence_correction') diff --git a/src/spectral_utilities.f90 b/src/spectral_utilities.f90 index dd753dc7c..5511de5cd 100644 --- a/src/spectral_utilities.f90 +++ b/src/spectral_utilities.f90 @@ -102,14 +102,6 @@ module spectral_utilities real(pReal) :: timeincOld end type tSolutionParams - type, public :: phaseFieldDataBin !< set of parameters defining a phase field - real(pReal) :: diffusion = 0.0_pReal, & !< thermal conductivity - mobility = 0.0_pReal, & !< thermal mobility - phaseField0 = 0.0_pReal !< homogeneous damage field starting condition - logical :: active = .false. - character(len=64) :: label = '' - end type phaseFieldDataBin - enum, bind(c) enumerator :: DERIVATIVE_CONTINUOUS_ID, & DERIVATIVE_CENTRAL_DIFF_ID, & From 4ee484b6e727579f3e4be413ea67aad9bde00b9b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 7 Mar 2019 07:12:35 +0100 Subject: [PATCH 3/4] correct scaling for fluctuations --- src/spectral_utilities.f90 | 244 ++++++++++++++++++------------------- 1 file changed, 118 insertions(+), 126 deletions(-) diff --git a/src/spectral_utilities.f90 b/src/spectral_utilities.f90 index 5511de5cd..532ec40b5 100644 --- a/src/spectral_utilities.f90 +++ b/src/spectral_utilities.f90 @@ -150,15 +150,9 @@ contains !> Initializes FFTW. !-------------------------------------------------------------------------------------------------- subroutine utilities_init() -#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 - use, intrinsic :: iso_fortran_env, only: & - compiler_version, & - compiler_options -#endif use IO, only: & IO_error, & IO_warning, & - IO_timeStamp, & IO_open_file use numerics, only: & spectral_derivative, & @@ -203,8 +197,6 @@ subroutine utilities_init() write(6,'(/,a)') ' <<<+- spectral_utilities init -+>>>' write(6,'(/,a)') ' Eisenlohr et al., International Journal of Plasticity, 46:37–53, 2013' write(6,'(a,/)') ' https://doi.org/10.1016/j.ijplas.2012.09.012' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() -#include "compilation_info.f90" !-------------------------------------------------------------------------------------------------- ! set debugging parameters @@ -1033,31 +1025,31 @@ end function utilities_calculateRate !> ensures that the average matches the aim !-------------------------------------------------------------------------------------------------- function utilities_forwardField(timeinc,field_lastInc,rate,aim) - use mesh, only: & - grid3, & - grid + use mesh, only: & + grid3, & + grid - implicit none - real(pReal), intent(in) :: & - timeinc !< timeinc of current step - real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: & - field_lastInc, & !< initial field - rate !< rate by which to forward - real(pReal), intent(in), optional, dimension(3,3) :: & - aim !< average field value aim - real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: & - utilities_forwardField - real(pReal), dimension(3,3) :: fieldDiff !< - aim - PetscErrorCode :: ierr + implicit none + real(pReal), intent(in) :: & + timeinc !< timeinc of current step + real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: & + field_lastInc, & !< initial field + rate !< rate by which to forward + real(pReal), intent(in), optional, dimension(3,3) :: & + aim !< average field value aim + real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: & + utilities_forwardField + real(pReal), dimension(3,3) :: fieldDiff !< - aim + PetscErrorCode :: ierr - utilities_forwardField = field_lastInc + rate*timeinc - if (present(aim)) then !< correct to match average - fieldDiff = sum(sum(sum(utilities_forwardField,dim=5),dim=4),dim=3)*wgt - call MPI_Allreduce(MPI_IN_PLACE,fieldDiff,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - fieldDiff = fieldDiff - aim - utilities_forwardField = utilities_forwardField - & - spread(spread(spread(fieldDiff,3,grid(1)),4,grid(2)),5,grid3) - endif + utilities_forwardField = field_lastInc + rate*timeinc + if (present(aim)) then !< correct to match average + fieldDiff = sum(sum(sum(utilities_forwardField,dim=5),dim=4),dim=3)*wgt + call MPI_Allreduce(MPI_IN_PLACE,fieldDiff,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + fieldDiff = fieldDiff - aim + utilities_forwardField = utilities_forwardField - & + spread(spread(spread(fieldDiff,3,grid(1)),4,grid(2)),5,grid3) + endif end function utilities_forwardField @@ -1068,51 +1060,50 @@ end function utilities_forwardField ! standard approach !-------------------------------------------------------------------------------------------------- pure function utilities_getFreqDerivative(k_s) - use math, only: & - PI - use mesh, only: & - geomSize, & - grid + use math, only: & + PI + use mesh, only: & + geomSize, & + grid - implicit none - integer(pInt), intent(in), dimension(3) :: k_s !< indices of frequency - complex(pReal), dimension(3) :: utilities_getFreqDerivative + implicit none + integer(pInt), intent(in), dimension(3) :: k_s !< indices of frequency + complex(pReal), dimension(3) :: utilities_getFreqDerivative - select case (spectral_derivative_ID) - case (DERIVATIVE_CONTINUOUS_ID) - utilities_getFreqDerivative = cmplx(0.0_pReal, 2.0_pReal*PI*real(k_s,pReal)/geomSize,pReal) + select case (spectral_derivative_ID) + case (DERIVATIVE_CONTINUOUS_ID) + utilities_getFreqDerivative = cmplx(0.0_pReal, 2.0_pReal*PI*real(k_s,pReal)/geomSize,pReal) - case (DERIVATIVE_CENTRAL_DIFF_ID) - utilities_getFreqDerivative = cmplx(0.0_pReal, sin(2.0_pReal*PI*real(k_s,pReal)/real(grid,pReal)), pReal)/ & - cmplx(2.0_pReal*geomSize/real(grid,pReal), 0.0_pReal, pReal) - - case (DERIVATIVE_FWBW_DIFF_ID) - utilities_getFreqDerivative(1) = & - cmplx(cos(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)) - 1.0_pReal, & - sin(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)), pReal)* & - cmplx(cos(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)) + 1.0_pReal, & - sin(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)), pReal)* & - cmplx(cos(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)) + 1.0_pReal, & - sin(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)), pReal)/ & - cmplx(4.0_pReal*geomSize(1)/real(grid(1),pReal), 0.0_pReal, pReal) - utilities_getFreqDerivative(2) = & - cmplx(cos(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)) + 1.0_pReal, & - sin(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)), pReal)* & - cmplx(cos(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)) - 1.0_pReal, & - sin(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)), pReal)* & - cmplx(cos(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)) + 1.0_pReal, & - sin(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)), pReal)/ & - cmplx(4.0_pReal*geomSize(2)/real(grid(2),pReal), 0.0_pReal, pReal) - utilities_getFreqDerivative(3) = & - cmplx(cos(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)) + 1.0_pReal, & - sin(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)), pReal)* & - cmplx(cos(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)) + 1.0_pReal, & - sin(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)), pReal)* & - cmplx(cos(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)) - 1.0_pReal, & - sin(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)), pReal)/ & - cmplx(4.0_pReal*geomSize(3)/real(grid(3),pReal), 0.0_pReal, pReal) - - end select + case (DERIVATIVE_CENTRAL_DIFF_ID) + utilities_getFreqDerivative = cmplx(0.0_pReal, sin(2.0_pReal*PI*real(k_s,pReal)/real(grid,pReal)), pReal)/ & + cmplx(2.0_pReal*geomSize/real(grid,pReal), 0.0_pReal, pReal) + + case (DERIVATIVE_FWBW_DIFF_ID) + utilities_getFreqDerivative(1) = & + cmplx(cos(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)) - 1.0_pReal, & + sin(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)), pReal)* & + cmplx(cos(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)) + 1.0_pReal, & + sin(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)), pReal)* & + cmplx(cos(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)) + 1.0_pReal, & + sin(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)), pReal)/ & + cmplx(4.0_pReal*geomSize(1)/real(grid(1),pReal), 0.0_pReal, pReal) + utilities_getFreqDerivative(2) = & + cmplx(cos(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)) + 1.0_pReal, & + sin(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)), pReal)* & + cmplx(cos(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)) - 1.0_pReal, & + sin(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)), pReal)* & + cmplx(cos(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)) + 1.0_pReal, & + sin(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)), pReal)/ & + cmplx(4.0_pReal*geomSize(2)/real(grid(2),pReal), 0.0_pReal, pReal) + utilities_getFreqDerivative(3) = & + cmplx(cos(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)) + 1.0_pReal, & + sin(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)), pReal)* & + cmplx(cos(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)) + 1.0_pReal, & + sin(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)), pReal)* & + cmplx(cos(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)) - 1.0_pReal, & + sin(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)), pReal)/ & + cmplx(4.0_pReal*geomSize(3)/real(grid(3),pReal), 0.0_pReal, pReal) + end select end function utilities_getFreqDerivative @@ -1123,59 +1114,60 @@ end function utilities_getFreqDerivative ! convolution !-------------------------------------------------------------------------------------------------- subroutine utilities_updateIPcoords(F) - use prec, only: & - cNeq - use IO, only: & - IO_error - use math, only: & - math_mul33x3 - use mesh, only: & - grid, & - grid3, & - grid3Offset, & - geomSize, & - mesh_ipCoordinates - implicit none - - real(pReal), dimension(3,3,grid(1),grid(2),grid3), intent(in) :: F - integer(pInt) :: i, j, k, m, ierr - real(pReal), dimension(3) :: step, offset_coords - real(pReal), dimension(3,3) :: Favg - -!-------------------------------------------------------------------------------------------------- -! integration in Fourier space - tensorField_real = 0.0_pReal - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F - call utilities_FFTtensorForward() - call utilities_fourierTensorDivergence() - - do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid1Red - if (any(cNeq(xi1st(1:3,i,j,k),cmplx(0.0_pReal,0.0_pReal)))) & - vectorField_fourier(1:3,i,j,k) = vectorField_fourier(1:3,i,j,k)/ & - sum(conjg(-xi1st(1:3,i,j,k))*xi1st(1:3,i,j,k)) - enddo; enddo; enddo - call fftw_mpi_execute_dft_c2r(planVectorBack,vectorField_fourier,vectorField_real) - -!-------------------------------------------------------------------------------------------------- -! average F - if (grid3Offset == 0_pInt) Favg = real(tensorField_fourier(1:3,1:3,1,1,1),pReal)*wgt - call MPI_Bcast(Favg,9,MPI_DOUBLE,0,PETSC_COMM_WORLD,ierr) - if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='update_IPcoords') - -!-------------------------------------------------------------------------------------------------- -! add average to fluctuation and put (0,0,0) on (0,0,0) - step = geomSize/real(grid, pReal) - if (grid3Offset == 0_pInt) offset_coords = vectorField_real(1:3,1,1,1) - call MPI_Bcast(offset_coords,3,MPI_DOUBLE,0,PETSC_COMM_WORLD,ierr) - if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='update_IPcoords') - offset_coords = math_mul33x3(Favg,step/2.0_pReal) - offset_coords - m = 1_pInt - do k = 1_pInt,grid3; do j = 1_pInt,grid(2); do i = 1_pInt,grid(1) - mesh_ipCoordinates(1:3,1,m) = vectorField_real(1:3,i,j,k) & - + offset_coords & - + math_mul33x3(Favg,step*real([i,j,k+grid3Offset]-1_pInt,pReal)) - m = m+1_pInt - enddo; enddo; enddo + use prec, only: & + cNeq + use IO, only: & + IO_error + use math, only: & + math_mul33x3 + use mesh, only: & + grid, & + grid3, & + grid3Offset, & + geomSize, & + mesh_ipCoordinates + implicit none + + real(pReal), dimension(3,3,grid(1),grid(2),grid3), intent(in) :: F + integer(pInt) :: i, j, k, m, ierr + real(pReal), dimension(3) :: step, offset_coords + real(pReal), dimension(3,3) :: Favg + + !-------------------------------------------------------------------------------------------------- + ! integration in Fourier space + tensorField_real = 0.0_pReal + tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F + call utilities_FFTtensorForward() + call utilities_fourierTensorDivergence() + + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid1Red + if (any(cNeq(xi1st(1:3,i,j,k),cmplx(0.0_pReal,0.0_pReal)))) & + vectorField_fourier(1:3,i,j,k) = vectorField_fourier(1:3,i,j,k)/ & + sum(conjg(-xi1st(1:3,i,j,k))*xi1st(1:3,i,j,k)) + enddo; enddo; enddo + call fftw_mpi_execute_dft_c2r(planVectorBack,vectorField_fourier,vectorField_real) + vectorField_real = vectorField_real * wgt + + !-------------------------------------------------------------------------------------------------- + ! average F + if (grid3Offset == 0_pInt) Favg = real(tensorField_fourier(1:3,1:3,1,1,1),pReal)*wgt + call MPI_Bcast(Favg,9,MPI_DOUBLE,0,PETSC_COMM_WORLD,ierr) + if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='update_IPcoords') + + !-------------------------------------------------------------------------------------------------- + ! add average to fluctuation and put (0,0,0) on (0,0,0) + step = geomSize/real(grid, pReal) + if (grid3Offset == 0_pInt) offset_coords = vectorField_real(1:3,1,1,1) + call MPI_Bcast(offset_coords,3,MPI_DOUBLE,0,PETSC_COMM_WORLD,ierr) + if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='update_IPcoords') + offset_coords = math_mul33x3(Favg,step/2.0_pReal) - offset_coords + m = 1_pInt + do k = 1_pInt,grid3; do j = 1_pInt,grid(2); do i = 1_pInt,grid(1) + mesh_ipCoordinates(1:3,1,m) = vectorField_real(1:3,i,j,k) & + + offset_coords & + + math_mul33x3(Favg,step*real([i,j,k+grid3Offset]-1_pInt,pReal)) + m = m+1_pInt + enddo; enddo; enddo end subroutine utilities_updateIPcoords From 4785f8fae468f4a0267a980857a4a72595614c08 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 7 Mar 2019 07:21:31 +0100 Subject: [PATCH 4/4] [skip ci] not needed anymore --- installation/patch/fwbw_derivative | 13 ------------- 1 file changed, 13 deletions(-) delete mode 100644 installation/patch/fwbw_derivative diff --git a/installation/patch/fwbw_derivative b/installation/patch/fwbw_derivative deleted file mode 100644 index 03d5be1c6..000000000 --- a/installation/patch/fwbw_derivative +++ /dev/null @@ -1,13 +0,0 @@ -diff --git a/code/numerics.f90 b/code/numerics.f90 -index 24bd190..c968c70 100644 ---- a/code/numerics.f90 -+++ b/code/numerics.f90 -@@ -110,7 +110,7 @@ module numerics - fftw_plan_mode = 'FFTW_PATIENT' !< reads the planing-rigor flag, see manual on www.fftw.org, Default FFTW_PATIENT: use patient planner flag - character(len=64), protected, public :: & - spectral_solver = 'basicpetsc' , & !< spectral solution method -- spectral_derivative = 'continuous' !< spectral spatial derivative method -+ spectral_derivative = 'fwbw_difference' !< spectral spatial derivative method - character(len=1024), protected, public :: & - petsc_defaultOptions = '-mech_snes_type ngmres & - &-damage_snes_type ngmres &