diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index 1b785ca3d..c9dea0166 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -106,8 +106,6 @@ program DAMASK_grid external :: & quit - class(tNode), pointer :: & - tmp type(tDict), pointer :: & config_load, & num_grid, & diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index 816cdd4b1..3e514b960 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -253,11 +253,13 @@ end function grid_damage_spectral_solution subroutine grid_damage_spectral_forward(cutBack) logical, intent(in) :: cutBack + integer :: i, j, k, ce DM :: dm_local PetscScalar, dimension(:,:,:), pointer :: phi_PETSc PetscErrorCode :: err_PETSc + if (cutBack) then phi_current = phi_lastInc phi_stagInc = phi_lastInc @@ -284,7 +286,7 @@ end subroutine grid_damage_spectral_forward !-------------------------------------------------------------------------------------------------- -!> @brief forms the spectral damage residual vector +!> @brief Construct the residual vector. !-------------------------------------------------------------------------------------------------- subroutine formResidual(in,x_scal,r,dummy,err_PETSc) @@ -297,48 +299,34 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc) X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: & r PetscObject :: dummy - PetscErrorCode :: err_PETSc + PetscErrorCode, intent(out) :: err_PETSc + integer :: i, j, k, ce + real(pReal), dimension(3,cells(1),cells(2),cells3) :: vectorField phi_current = x_scal !-------------------------------------------------------------------------------------------------- ! evaluate polarization field - scalarField_real = 0.0_pReal - scalarField_real(1:cells(1),1:cells(2),1:cells3) = phi_current - call utilities_FFTscalarForward - call utilities_fourierScalarGradient !< calculate gradient of damage field - call utilities_FFTvectorBackward + vectorField = utilities_ScalarGradient(phi_current) ce = 0 do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) ce = ce + 1 - vectorField_real(1:3,i,j,k) = matmul(homogenization_K_phi(ce) - K_ref, vectorField_real(1:3,i,j,k)) + vectorField(1:3,i,j,k) = matmul(homogenization_K_phi(ce) - K_ref, vectorField(1:3,i,j,k)) end do; end do; end do - call utilities_FFTvectorForward - call utilities_fourierVectorDivergence !< calculate damage divergence in fourier field - call utilities_FFTscalarBackward + r = utilities_VectorDivergence(vectorField) ce = 0 do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) ce = ce + 1 - scalarField_real(i,j,k) = params%Delta_t*(scalarField_real(i,j,k) + homogenization_f_phi(phi_current(i,j,k),ce)) & - + homogenization_mu_phi(ce)*(phi_lastInc(i,j,k) - phi_current(i,j,k)) & - + mu_ref*phi_current(i,j,k) + r(i,j,k) = params%Delta_t*(r(i,j,k) + homogenization_f_phi(phi_current(i,j,k),ce)) & + + homogenization_mu_phi(ce)*(phi_lastInc(i,j,k) - phi_current(i,j,k)) & + + mu_ref*phi_current(i,j,k) end do; end do; end do -!-------------------------------------------------------------------------------------------------- -! convolution of damage field with green operator - call utilities_FFTscalarForward - call utilities_fourierGreenConvolution(K_ref, mu_ref, params%Delta_t) - call utilities_FFTscalarBackward - - where(scalarField_real(1:cells(1),1:cells(2),1:cells3) > phi_lastInc) & - scalarField_real(1:cells(1),1:cells(2),1:cells3) = phi_lastInc - where(scalarField_real(1:cells(1),1:cells(2),1:cells3) < num%residualStiffness) & - scalarField_real(1:cells(1),1:cells(2),1:cells3) = num%residualStiffness - !-------------------------------------------------------------------------------------------------- ! constructing residual - r = scalarField_real(1:cells(1),1:cells(2),1:cells3) - phi_current + r = max(min(utilities_GreenConvolution(r, K_ref, mu_ref, params%Delta_t),phi_lastInc),num%residualStiffness) & + - phi_current err_PETSc = 0 end subroutine formResidual diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index bead280a7..c570e94e7 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -491,7 +491,7 @@ end subroutine converged !-------------------------------------------------------------------------------------------------- -!> @brief forms the residual vector +!> @brief Construct the residual vector. !-------------------------------------------------------------------------------------------------- subroutine formResidual(in, F, & r, dummy, err_PETSc) @@ -501,15 +501,17 @@ subroutine formResidual(in, F, & intent(in) :: F !< deformation gradient field PetscScalar, dimension(3,3,X_RANGE,Y_RANGE,Z_RANGE), & intent(out) :: r !< residuum field + PetscObject :: dummy + PetscErrorCode :: err_PETSc + real(pReal), dimension(3,3) :: & deltaF_aim PetscInt :: & PETScIter, & nfuncs - PetscObject :: dummy - PetscErrorCode :: err_PETSc integer(MPI_INTEGER_KIND) :: err_MPI + call SNESGetNumberFunctionEvals(SNES_mechanical,nfuncs,err_PETSc) CHKERRQ(err_PETSc) call SNESGetIterationNumber(SNES_mechanical,PETScIter,err_PETSc) @@ -517,8 +519,6 @@ subroutine formResidual(in, F, & if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment -!-------------------------------------------------------------------------------------------------- -! begin of new iteration newIteration: if (totalIter <= PETScIter) then totalIter = totalIter + 1 print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax @@ -529,32 +529,20 @@ subroutine formResidual(in, F, & flush(IO_STDOUT) end if newIteration -!-------------------------------------------------------------------------------------------------- -! evaluate constitutive response - call utilities_constitutiveResponse(r, & ! residuum gets field of first PK stress (to save memory) - P_av,C_volAvg,C_minMaxAvg, & - F,params%Delta_t,params%rotation_BC) - call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + associate (P => r) + call utilities_constitutiveResponse(P, & + P_av,C_volAvg,C_minMaxAvg, & + F,params%Delta_t,params%rotation_BC) + call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) + if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + err_div = utilities_divergenceRMS(P) + end associate -!-------------------------------------------------------------------------------------------------- -! stress BC handling deltaF_aim = math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc F_aim = F_aim - deltaF_aim err_BC = maxval(abs(merge(.0_pReal,P_av - P_aim,params%stress_mask))) -!-------------------------------------------------------------------------------------------------- -! updated deformation gradient using fix point algorithm of basic scheme - tensorField_real = 0.0_pReal - tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = r ! store fPK field for subsequent FFT forward transform - call utilities_FFTtensorForward ! FFT forward of global "tensorField_real" - err_div = utilities_divergenceRMS() ! divRMS of tensorField_fourier for later use - call utilities_fourierGammaConvolution(params%rotation_BC%rotate(deltaF_aim,active=.true.)) ! convolution of Gamma and tensorField_fourier - call utilities_FFTtensorBackward ! FFT backward of global tensorField_fourier - -!-------------------------------------------------------------------------------------------------- -! constructing residual - r = tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too + r = utilities_GammaConvolution(r,params%rotation_BC%rotate(deltaF_aim,active=.true.)) end subroutine formResidual diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 2b4ea364a..46ff04adb 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -551,7 +551,7 @@ end subroutine converged !-------------------------------------------------------------------------------------------------- -!> @brief forms the residual vector +!> @brief Construct the residual vector. !-------------------------------------------------------------------------------------------------- subroutine formResidual(in, FandF_tau, & r, dummy,err_PETSc) @@ -561,6 +561,9 @@ subroutine formResidual(in, FandF_tau, & target, intent(in) :: FandF_tau PetscScalar, dimension(3,3,2,X_RANGE,Y_RANGE,Z_RANGE),& target, intent(out) :: r !< residuum field + PetscObject :: dummy + PetscErrorCode :: err_PETSc + PetscScalar, pointer, dimension(:,:,:,:,:) :: & F, & F_tau, & @@ -569,13 +572,10 @@ subroutine formResidual(in, FandF_tau, & PetscInt :: & PETScIter, & nfuncs - PetscObject :: dummy - PetscErrorCode :: err_PETSc integer(MPI_INTEGER_KIND) :: err_MPI integer :: & i, j, k, e -!--------------------------------------------------------------------------------------------------- F => FandF_tau(1:3,1:3,1,& XG_RANGE,YG_RANGE,ZG_RANGE) @@ -597,8 +597,6 @@ subroutine formResidual(in, FandF_tau, & if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment -!-------------------------------------------------------------------------------------------------- -! begin of new iteration newIteration: if (totalIter <= PETScIter) then totalIter = totalIter + 1 print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax @@ -609,63 +607,53 @@ subroutine formResidual(in, FandF_tau, & flush(IO_STDOUT) end if newIteration -!-------------------------------------------------------------------------------------------------- -! - tensorField_real = 0.0_pReal do k = 1, cells3; do j = 1, cells(2); do i = 1, cells(1) - tensorField_real(1:3,1:3,i,j,k) = & + r_F_tau(1:3,1:3,i,j,k) = & num%beta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -& num%alpha*matmul(F(1:3,1:3,i,j,k), & math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3)) end do; end do; end do + r_F_tau = num%beta*F & + - utilities_GammaConvolution(r_F_tau,params%rotation_BC%rotate(num%beta*F_aim,active=.true.)) -!-------------------------------------------------------------------------------------------------- -! doing convolution in Fourier space - call utilities_FFTtensorForward - call utilities_fourierGammaConvolution(params%rotation_BC%rotate(num%beta*F_aim,active=.true.)) - call utilities_FFTtensorBackward + err_curl = utilities_curlRMS(F) -!-------------------------------------------------------------------------------------------------- -! constructing residual - r_F_tau = num%beta*F - tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) +#ifdef __GFORTRAN__ + call utilities_constitutiveResponse(r_F, & +#else + associate (P => r_F) + call utilities_constitutiveResponse(P, & +#endif + P_av,C_volAvg,C_minMaxAvg, & + F - r_F_tau/num%beta,params%Delta_t,params%rotation_BC) + call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) +#ifdef __GFORTRAN__ + err_div = utilities_divergenceRMS(r_F) +#else + err_div = utilities_divergenceRMS(P) +#endif + e = 0 + do k = 1, cells3; do j = 1, cells(2); do i = 1, cells(1) + e = e + 1 + r_F(1:3,1:3,i,j,k) = & + math_mul3333xx33(math_invSym3333(homogenization_dPdF(1:3,1:3,1:3,1:3,e) + C_scale), & +#ifdef __GFORTRAN__ + r_F(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), & +#else + P(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), & +#endif + math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))) & + + r_F_tau(1:3,1:3,i,j,k) + end do; end do; end do +#ifndef __GFORTRAN__ + end associate +#endif -!-------------------------------------------------------------------------------------------------- -! evaluate constitutive response - call utilities_constitutiveResponse(r_F, & ! "residuum" gets field of first PK stress (to save memory) - P_av,C_volAvg,C_minMaxAvg, & - F - r_F_tau/num%beta,params%Delta_t,params%rotation_BC) - call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) - -!-------------------------------------------------------------------------------------------------- -! stress BC handling F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc err_BC = maxval(abs(merge(math_mul3333xx33(C_scale,F_aim-params%rotation_BC%rotate(F_av)), & P_av-P_aim, & params%stress_mask))) -! calculate divergence - tensorField_real = 0.0_pReal - tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = r_F !< stress field in disguise - call utilities_FFTtensorForward - err_div = utilities_divergenceRMS() !< root mean squared error in divergence of stress -!-------------------------------------------------------------------------------------------------- -! constructing residual - e = 0 - do k = 1, cells3; do j = 1, cells(2); do i = 1, cells(1) - e = e + 1 - r_F(1:3,1:3,i,j,k) = & - math_mul3333xx33(math_invSym3333(homogenization_dPdF(1:3,1:3,1:3,1:3,e) + C_scale), & - r_F(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), & - math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))) & - + r_F_tau(1:3,1:3,i,j,k) - end do; end do; end do - -!-------------------------------------------------------------------------------------------------- -! calculating curl - tensorField_real = 0.0_pReal - tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = F - call utilities_FFTtensorForward - err_curl = utilities_curlRMS() end subroutine formResidual diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 29ae07769..b1dde568f 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -242,11 +242,13 @@ end function grid_thermal_spectral_solution subroutine grid_thermal_spectral_forward(cutBack) logical, intent(in) :: cutBack + integer :: i, j, k, ce DM :: dm_local PetscScalar, dimension(:,:,:), pointer :: T_PETSc PetscErrorCode :: err_PETSc + if (cutBack) then T_current = T_lastInc T_stagInc = T_lastInc @@ -307,7 +309,7 @@ end subroutine grid_thermal_spectral_restartWrite !-------------------------------------------------------------------------------------------------- -!> @brief forms the spectral thermal residual vector +!> @brief Construct the residual vector. !-------------------------------------------------------------------------------------------------- subroutine formResidual(in,x_scal,r,dummy,err_PETSc) @@ -320,42 +322,34 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc) X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: & r PetscObject :: dummy - PetscErrorCode :: err_PETSc + PetscErrorCode, intent(out) :: err_PETSc + integer :: i, j, k, ce + real(pReal), dimension(3,cells(1),cells(2),cells3) :: vectorField + T_current = x_scal !-------------------------------------------------------------------------------------------------- ! evaluate polarization field - scalarField_real = 0.0_pReal - scalarField_real(1:cells(1),1:cells(2),1:cells3) = T_current - call utilities_FFTscalarForward - call utilities_fourierScalarGradient !< calculate gradient of temperature field - call utilities_FFTvectorBackward + vectorField = utilities_ScalarGradient(T_current) ce = 0 do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) ce = ce + 1 - vectorField_real(1:3,i,j,k) = matmul(homogenization_K_T(ce) - K_ref, vectorField_real(1:3,i,j,k)) + vectorField(1:3,i,j,k) = matmul(homogenization_K_T(ce) - K_ref, vectorField(1:3,i,j,k)) end do; end do; end do - call utilities_FFTvectorForward - call utilities_fourierVectorDivergence !< calculate temperature divergence in fourier field - call utilities_FFTscalarBackward + r = utilities_VectorDivergence(vectorField) ce = 0 do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) ce = ce + 1 - scalarField_real(i,j,k) = params%Delta_t*(scalarField_real(i,j,k) + homogenization_f_T(ce)) & - + homogenization_mu_T(ce) * (T_lastInc(i,j,k) - T_current(i,j,k)) & - + mu_ref*T_current(i,j,k) + r(i,j,k) = params%Delta_t*(r(i,j,k) + homogenization_f_T(ce)) & + + homogenization_mu_T(ce) * (T_lastInc(i,j,k) - T_current(i,j,k)) & + + mu_ref*T_current(i,j,k) end do; end do; end do -!-------------------------------------------------------------------------------------------------- -! convolution of temperature field with green operator - call utilities_FFTscalarForward - call utilities_fourierGreenConvolution(K_ref, mu_ref, params%Delta_t) - call utilities_FFTscalarBackward - !-------------------------------------------------------------------------------------------------- ! constructing residual - r = T_current - scalarField_real(1:cells(1),1:cells(2),1:cells3) + r = T_current & + - utilities_GreenConvolution(r, K_ref, mu_ref, params%Delta_t) err_PETSc = 0 end subroutine formResidual diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 6cb7edd30..54dbf8c43 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -42,16 +42,16 @@ module spectral_utilities !-------------------------------------------------------------------------------------------------- ! variables storing information for spectral method and FFTW - real(C_DOUBLE), public, dimension(:,:,:,:,:), pointer :: tensorField_real !< tensor field in real space - real(C_DOUBLE), public, dimension(:,:,:,:), pointer :: vectorField_real !< vector field in real space - real(C_DOUBLE), public, dimension(:,:,:), pointer :: scalarField_real !< scalar field in real space - complex(C_DOUBLE_COMPLEX), dimension(:,:,:,:,:), pointer :: tensorField_fourier !< tensor field in Fourier space - complex(C_DOUBLE_COMPLEX), dimension(:,:,:,:), pointer :: vectorField_fourier !< vector field in Fourier space - complex(C_DOUBLE_COMPLEX), dimension(:,:,:), pointer :: scalarField_fourier !< scalar field in Fourier space - complex(pReal), dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat !< gamma operator (field) for spectral method - complex(pReal), dimension(:,:,:,:), allocatable :: xi1st !< wave vector field for first derivatives - complex(pReal), dimension(:,:,:,:), allocatable :: xi2nd !< wave vector field for second derivatives - real(pReal), dimension(3,3,3,3) :: C_ref !< mechanic reference stiffness + real(C_DOUBLE), dimension(:,:,:,:,:), pointer :: tensorField_real !< tensor field in real space + real(C_DOUBLE), dimension(:,:,:,:), pointer :: vectorField_real !< vector field in real space + real(C_DOUBLE), dimension(:,:,:), pointer :: scalarField_real !< scalar field in real space + complex(C_DOUBLE_COMPLEX), dimension(:,:,:,:,:), pointer :: tensorField_fourier !< tensor field in Fourier space + complex(C_DOUBLE_COMPLEX), dimension(:,:,:,:), pointer :: vectorField_fourier !< vector field in Fourier space + complex(C_DOUBLE_COMPLEX), dimension(:,:,:), pointer :: scalarField_fourier !< scalar field in Fourier space + complex(pReal), dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat !< gamma operator (field) for spectral method + complex(pReal), dimension(:,:,:,:), allocatable :: xi1st !< wave vector field for first derivatives + complex(pReal), dimension(:,:,:,:), allocatable :: xi2nd !< wave vector field for second derivatives + real(pReal), dimension(3,3,3,3) :: C_ref !< mechanic reference stiffness !-------------------------------------------------------------------------------------------------- @@ -116,18 +116,12 @@ module spectral_utilities public :: & spectral_utilities_init, & utilities_updateGamma, & - utilities_FFTtensorForward, & - utilities_FFTtensorBackward, & - utilities_FFTvectorForward, & - utilities_FFTvectorBackward, & - utilities_FFTscalarForward, & - utilities_FFTscalarBackward, & - utilities_fourierGammaConvolution, & - utilities_fourierGreenConvolution, & + utilities_GammaConvolution, & + utilities_GreenConvolution, & utilities_divergenceRMS, & utilities_curlRMS, & - utilities_fourierScalarGradient, & - utilities_fourierVectorDivergence, & + utilities_ScalarGradient, & + utilities_VectorDivergence, & utilities_maskedCompliance, & utilities_constitutiveResponse, & utilities_calculateRate, & @@ -306,8 +300,8 @@ subroutine spectral_utilities_init() !-------------------------------------------------------------------------------------------------- ! allocation - allocate (xi1st (3,cells1Red,cells(3),cells2),source = cmplx(0.0_pReal,0.0_pReal,pReal)) ! frequencies for first derivatives, only half the size for first dimension - allocate (xi2nd (3,cells1Red,cells(3),cells2),source = cmplx(0.0_pReal,0.0_pReal,pReal)) ! frequencies for second derivatives, only half the size for first dimension + allocate (xi1st (3,cells1Red,cells(3),cells2),source = cmplx(0.0_pReal,0.0_pReal,pReal)) ! frequencies for first derivatives, only half the size for first dimension + allocate (xi2nd (3,cells1Red,cells(3),cells2),source = cmplx(0.0_pReal,0.0_pReal,pReal)) ! frequencies for second derivatives, only half the size for first dimension !-------------------------------------------------------------------------------------------------- ! tensor MPI fftw plans @@ -385,6 +379,7 @@ end subroutine spectral_utilities_init subroutine utilities_updateGamma(C) real(pReal), intent(in), dimension(3,3,3,3) :: C !< input stiffness to store as reference stiffness + complex(pReal), dimension(3,3) :: temp33_cmplx, xiDyad_cmplx real(pReal), dimension(6,6) :: A, A_inv integer :: & @@ -392,7 +387,8 @@ subroutine utilities_updateGamma(C) l, m, n, o logical :: err - C_ref = C + + C_ref = C/wgt if (.not. num%memory_efficient) then gamma_hat = cmplx(0.0_pReal,0.0_pReal,pReal) ! for the singular point and any non invertible A @@ -434,68 +430,6 @@ subroutine utilities_updateGamma(C) end subroutine utilities_updateGamma -!-------------------------------------------------------------------------------------------------- -!> @brief forward FFT of data in field_real to field_fourier -!> @details Does an unweighted FFT transform from real to complex. Extra padding entries are set -! to 0.0 -!-------------------------------------------------------------------------------------------------- -subroutine utilities_FFTtensorForward() - - tensorField_real(1:3,1:3,cells(1)+1:cells1Red*2,:,:) = 0.0_pReal - call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier) - -end subroutine utilities_FFTtensorForward - - -!-------------------------------------------------------------------------------------------------- -!> @brief backward FFT of data in field_fourier to field_real -!> @details Does an weighted inverse FFT transform from complex to real -!-------------------------------------------------------------------------------------------------- -subroutine utilities_FFTtensorBackward() - - call fftw_mpi_execute_dft_c2r(planTensorBack,tensorField_fourier,tensorField_real) - tensorField_real = tensorField_real * wgt ! normalize the result by number of elements - -end subroutine utilities_FFTtensorBackward - -!-------------------------------------------------------------------------------------------------- -!> @brief forward FFT of data in scalarField_real to scalarField_fourier -!> @details Does an unweighted FFT transform from real to complex. Extra padding entries are set -! to 0.0 -!-------------------------------------------------------------------------------------------------- -subroutine utilities_FFTscalarForward() - - scalarField_real(cells(1)+1:cells1Red*2,:,:) = 0.0_pReal - call fftw_mpi_execute_dft_r2c(planScalarForth,scalarField_real,scalarField_fourier) - -end subroutine utilities_FFTscalarForward - - -!-------------------------------------------------------------------------------------------------- -!> @brief backward FFT of data in scalarField_fourier to scalarField_real -!> @details Does an weighted inverse FFT transform from complex to real -!-------------------------------------------------------------------------------------------------- -subroutine utilities_FFTscalarBackward() - - call fftw_mpi_execute_dft_c2r(planScalarBack,scalarField_fourier,scalarField_real) - scalarField_real = scalarField_real * wgt ! normalize the result by number of elements - -end subroutine utilities_FFTscalarBackward - - -!-------------------------------------------------------------------------------------------------- -!> @brief forward FFT of data in field_real to field_fourier with highest freqs. removed -!> @details Does an unweighted FFT transform from real to complex. Extra padding entries are set -! to 0.0 -!-------------------------------------------------------------------------------------------------- -subroutine utilities_FFTvectorForward() - - vectorField_real(1:3,cells(1)+1:cells1Red*2,:,:) = 0.0_pReal - call fftw_mpi_execute_dft_r2c(planVectorForth,vectorField_real,vectorField_fourier) - -end subroutine utilities_FFTvectorForward - - !-------------------------------------------------------------------------------------------------- !> @brief backward FFT of data in field_fourier to field_real !> @details Does an weighted inverse FFT transform from complex to real @@ -511,12 +445,14 @@ end subroutine utilities_FFTvectorBackward !-------------------------------------------------------------------------------------------------- !> @brief doing convolution gamma_hat * field_real, ensuring that average value = fieldAim !-------------------------------------------------------------------------------------------------- -subroutine utilities_fourierGammaConvolution(fieldAim) +function utilities_GammaConvolution(field, fieldAim) result(gammaField) + real(pReal), intent(in), dimension(3,3,cells(1),cells(2),cells3) :: field real(pReal), intent(in), dimension(3,3) :: fieldAim !< desired average value of the field after convolution - complex(pReal), dimension(3,3) :: temp33_cmplx, xiDyad_cmplx - real(pReal), dimension(6,6) :: A, A_inv + real(pReal), dimension(3,3,cells(1),cells(2),cells3) :: gammaField + complex(pReal), dimension(3,3) :: temp33_cmplx, xiDyad_cmplx + real(pReal), dimension(6,6) :: A, A_inv integer :: & i, j, k, & l, m, n, o @@ -526,8 +462,10 @@ subroutine utilities_fourierGammaConvolution(fieldAim) print'(/,1x,a)', '... doing gamma convolution ...............................................' flush(IO_STDOUT) -!-------------------------------------------------------------------------------------------------- -! do the actual spectral method calculation (mechanical equilibrium) + tensorField_real(1:3,1:3,cells(1)+1:cells1Red*2,1:cells(2),1:cells3) = 0.0_pReal + tensorField_real(1:3,1:3,1:cells(1), 1:cells(2),1:cells3) = field + call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier) + memoryEfficient: if (num%memory_efficient) then !$OMP PARALLEL DO PRIVATE(l,m,n,o,temp33_cmplx,xiDyad_cmplx,A,A_inv,err,gamma_hat) do j = 1, cells2; do k = 1, cells(3); do i = 1, cells1Red @@ -586,39 +524,53 @@ subroutine utilities_fourierGammaConvolution(fieldAim) !$OMP END PARALLEL DO end if memoryEfficient - if (cells3Offset == 0) tensorField_fourier(1:3,1:3,1,1,1) = cmplx(fieldAim/wgt,0.0_pReal,pReal) + if (cells3Offset == 0) tensorField_fourier(1:3,1:3,1,1,1) = cmplx(fieldAim,0.0_pReal,pReal) -end subroutine utilities_fourierGammaConvolution + call fftw_mpi_execute_dft_c2r(planTensorBack,tensorField_fourier,tensorField_real) + gammaField = tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) + +end function utilities_GammaConvolution !-------------------------------------------------------------------------------------------------- -!> @brief doing convolution DamageGreenOp_hat * field_real +!> @brief Convolution of Greens' operator for damage/thermal. !-------------------------------------------------------------------------------------------------- -subroutine utilities_fourierGreenConvolution(D_ref, mu_ref, Delta_t) +function utilities_GreenConvolution(field, D_ref, mu_ref, Delta_t) result(greenField) + real(pReal), intent(in), dimension(cells(1),cells(2),cells3) :: field real(pReal), dimension(3,3), intent(in) :: D_ref real(pReal), intent(in) :: mu_ref, Delta_t + real(pReal), dimension(cells(1),cells(2),cells3) :: greenField + complex(pReal) :: GreenOp_hat integer :: i, j, k -!-------------------------------------------------------------------------------------------------- -! do the actual spectral method calculation + + scalarField_real(cells(1)+1:cells1Red*2,1:cells(2),1:cells3) = 0.0_pReal + scalarField_real(1:cells(1), 1:cells(2),1:cells3) = field + call fftw_mpi_execute_dft_r2c(planScalarForth,scalarField_real,scalarField_fourier) + !$OMP PARALLEL DO PRIVATE(GreenOp_hat) do j = 1, cells2; do k = 1, cells(3); do i = 1, cells1Red - GreenOp_hat = cmplx(1.0_pReal,0.0_pReal,pReal) & + GreenOp_hat = cmplx(wgt,0.0_pReal,pReal) & / (cmplx(mu_ref,0.0_pReal,pReal) + cmplx(Delta_t,0.0_pReal,pReal) & * sum(conjg(xi1st(1:3,i,k,j))* matmul(cmplx(D_ref,0.0_pReal,pReal),xi1st(1:3,i,k,j)))) scalarField_fourier(i,k,j) = scalarField_fourier(i,k,j)*GreenOp_hat end do; end do; end do !$OMP END PARALLEL DO -end subroutine utilities_fourierGreenConvolution + call fftw_mpi_execute_dft_c2r(planScalarBack,scalarField_fourier,scalarField_real) + greenField = scalarField_real(1:cells(1),1:cells(2),1:cells3) + +end function utilities_GreenConvolution !-------------------------------------------------------------------------------------------------- -!> @brief calculate root mean square of divergence of field_fourier +!> @brief Calculate root mean square of divergence. !-------------------------------------------------------------------------------------------------- -real(pReal) function utilities_divergenceRMS() +real(pReal) function utilities_divergenceRMS(tensorField) + + real(pReal), dimension(3,3,cells(1),cells(2),cells3), intent(in) :: tensorField integer :: i, j, k integer(MPI_INTEGER_KIND) :: err_MPI @@ -628,6 +580,10 @@ real(pReal) function utilities_divergenceRMS() print'(/,1x,a)', '... calculating divergence ................................................' flush(IO_STDOUT) + tensorField_real(1:3,1:3,cells(1)+1:cells1Red*2,1:cells(2),1:cells3) = 0.0_pReal + tensorField_real(1:3,1:3,1:cells(1), 1:cells(2),1:cells3) = tensorField + call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier) + rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal,pReal) !-------------------------------------------------------------------------------------------------- @@ -660,9 +616,11 @@ end function utilities_divergenceRMS !-------------------------------------------------------------------------------------------------- -!> @brief calculate max of curl of field_fourier +!> @brief Calculate root mean square of curl. !-------------------------------------------------------------------------------------------------- -real(pReal) function utilities_curlRMS() +real(pReal) function utilities_curlRMS(tensorField) + + real(pReal), dimension(3,3,cells(1),cells(2),cells3), intent(in) :: tensorField integer :: i, j, k, l integer(MPI_INTEGER_KIND) :: err_MPI @@ -673,6 +631,10 @@ real(pReal) function utilities_curlRMS() print'(/,1x,a)', '... calculating curl ......................................................' flush(IO_STDOUT) + tensorField_real(1:3,1:3,cells(1)+1:cells1Red*2,1:cells(2),1:cells3) = 0.0_pReal + tensorField_real(1:3,1:3,1:cells(1), 1:cells(2),1:cells3) = tensorField + call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier) + rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal,pReal) !-------------------------------------------------------------------------------------------------- @@ -723,7 +685,7 @@ end function utilities_curlRMS !-------------------------------------------------------------------------------------------------- -!> @brief calculates mask compliance tensor used to adjust F to fullfill stress BC +!> @brief Calculate masked compliance tensor used to adjust F to fullfill stress BC. !-------------------------------------------------------------------------------------------------- function utilities_maskedCompliance(rot_BC,mask_stress,C) @@ -793,29 +755,46 @@ end function utilities_maskedCompliance !-------------------------------------------------------------------------------------------------- -!> @brief calculate scalar gradient in fourier field +!> @brief Calculate gradient of scalar field. !-------------------------------------------------------------------------------------------------- -subroutine utilities_fourierScalarGradient() +function utilities_ScalarGradient(field) result(grad) + + real(pReal), intent(in), dimension( cells(1),cells(2),cells3) :: field + real(pReal), dimension(3,cells(1),cells(2),cells3) :: grad integer :: i, j, k + scalarField_real(cells(1)+1:cells1Red*2,1:cells(2),1:cells3) = 0.0_pReal + scalarField_real(1:cells(1), 1:cells(2),1:cells3) = field + call fftw_mpi_execute_dft_r2c(planScalarForth,scalarField_real,scalarField_fourier) do j = 1, cells2; do k = 1, cells(3); do i = 1,cells1Red vectorField_fourier(1:3,i,k,j) = scalarField_fourier(i,k,j)*xi1st(1:3,i,k,j) ! ToDo: no -conjg? end do; end do; end do + call fftw_mpi_execute_dft_c2r(planVectorBack,vectorField_fourier,vectorField_real) + grad = vectorField_real(1:3,1:cells(1),1:cells(2),1:cells3)*wgt -end subroutine utilities_fourierScalarGradient +end function utilities_ScalarGradient !-------------------------------------------------------------------------------------------------- -!> @brief calculate vector divergence in fourier field +!> @brief Calculate divergence of vector field. !-------------------------------------------------------------------------------------------------- -subroutine utilities_fourierVectorDivergence() +function utilities_VectorDivergence(field) result(div) + real(pReal), intent(in), dimension(3,cells(1),cells(2),cells3) :: field + real(pReal), dimension( cells(1),cells(2),cells3) :: div + + + vectorField_real(1:3,cells(1)+1:cells1Red*2,1:cells(2),1:cells3) = 0.0_pReal + vectorField_real(1:3,1:cells(1), 1:cells(2),1:cells3) = field + call fftw_mpi_execute_dft_r2c(planVectorForth,vectorField_real,vectorField_fourier) scalarField_fourier(1:cells1Red,1:cells(3),1:cells2) = sum(vectorField_fourier(1:3,1:cells1Red,1:cells(3),1:cells2) & *conjg(-xi1st),1) + call fftw_mpi_execute_dft_c2r(planScalarBack,scalarField_fourier,scalarField_real) + div = scalarField_real(1:cells(1),1:cells(2),1:cells3)*wgt -end subroutine utilities_fourierVectorDivergence +end function utilities_VectorDivergence !-------------------------------------------------------------------------------------------------- @@ -1046,11 +1025,19 @@ subroutine utilities_updateCoords(F) step = geomSize/real(cells, pReal) + + tensorField_real(1:3,1:3,1:cells(1), 1:cells(2),1:cells3) = F + tensorField_real(1:3,1:3,cells(1)+1:cells1Red*2,1:cells(2),1:cells3) = 0.0_pReal + call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier) + + !-------------------------------------------------------------------------------------------------- + ! average F + if (cells3Offset == 0) Favg = tensorField_fourier(1:3,1:3,1,1,1)%re*wgt + call MPI_Bcast(Favg,9_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) + if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + !-------------------------------------------------------------------------------------------------- ! integration in Fourier space to get fluctuations of cell center discplacements - tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = F - call utilities_FFTtensorForward() - !$OMP PARALLEL DO do j = 1, cells2; do k = 1, cells(3); do i = 1, cells1Red if (any([i,j+cells2Offset,k] /= 1)) then @@ -1062,13 +1049,8 @@ subroutine utilities_updateCoords(F) end do; end do; end do !$OMP END PARALLEL DO - call utilities_FFTvectorBackward() - - !-------------------------------------------------------------------------------------------------- - ! average F - if (cells3Offset == 0) Favg = tensorField_fourier(1:3,1:3,1,1,1)%re*wgt - call MPI_Bcast(Favg,9_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + call fftw_mpi_execute_dft_c2r(planVectorBack,vectorField_fourier,vectorField_real) + vectorField_real = vectorField_real * wgt ! normalize the result by number of elements !-------------------------------------------------------------------------------------------------- ! pad cell center fluctuations along z-direction (needed when running MPI simulation) @@ -1136,7 +1118,7 @@ subroutine utilities_saveReferenceStiffness() open(newunit=fileUnit, file=getSolverJobName()//'.C_ref',& status='replace',access='stream',action='write',iostat=ierr) if (ierr /=0) call IO_error(100,ext_msg='could not open file '//getSolverJobName()//'.C_ref') - write(fileUnit) C_ref + write(fileUnit) C_ref*wgt close(fileUnit) end if @@ -1156,38 +1138,38 @@ subroutine selfTest() call random_number(tensorField_real) tensorField_real(1:3,1:3,cells(1)+1:cells1Red*2,:,:) = 0.0_pReal tensorField_real_ = tensorField_real - call utilities_FFTtensorForward() + call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier) if (worldsize==1) then if (any(dNeq(sum(sum(sum(tensorField_real_,dim=5),dim=4),dim=3)/tensorField_fourier(:,:,1,1,1)%re,1.0_pReal,1.0e-12_pReal))) & error stop 'tensorField avg' endif - call utilities_FFTtensorBackward() + call fftw_mpi_execute_dft_c2r(planTensorBack,tensorField_fourier,tensorField_real) tensorField_real(1:3,1:3,cells(1)+1:cells1Red*2,:,:) = 0.0_pReal - if (maxval(abs(tensorField_real_ - tensorField_real))>5.0e-15_pReal) error stop 'tensorField' + if (maxval(abs(tensorField_real_ - tensorField_real*wgt))>5.0e-15_pReal) error stop 'tensorField' call random_number(vectorField_real) vectorField_real(1:3,cells(1)+1:cells1Red*2,:,:) = 0.0_pReal vectorField_real_ = vectorField_real - call utilities_FFTvectorForward() + call fftw_mpi_execute_dft_r2c(planVectorForth,vectorField_real,vectorField_fourier) if (worldsize==1) then if (any(dNeq(sum(sum(sum(vectorField_real_,dim=4),dim=3),dim=2)/vectorField_fourier(:,1,1,1)%re,1.0_pReal,1.0e-12_pReal))) & error stop 'vector avg' endif - call utilities_FFTvectorBackward() + call fftw_mpi_execute_dft_c2r(planVectorBack,vectorField_fourier,vectorField_real) vectorField_real(1:3,cells(1)+1:cells1Red*2,:,:) = 0.0_pReal - if (maxval(abs(vectorField_real_ - vectorField_real))>5.0e-15_pReal) error stop 'vectorField' + if (maxval(abs(vectorField_real_ - vectorField_real*wgt))>5.0e-15_pReal) error stop 'vectorField' call random_number(scalarField_real) scalarField_real(cells(1)+1:cells1Red*2,:,:) = 0.0_pReal scalarField_real_ = scalarField_real - call utilities_FFTscalarForward() + call fftw_mpi_execute_dft_r2c(planScalarForth,scalarField_real,scalarField_fourier) if (worldsize==1) then if (dNeq(sum(sum(sum(scalarField_real_,dim=3),dim=2),dim=1)/scalarField_fourier(1,1,1)%re,1.0_pReal,1.0e-12_pReal)) & error stop 'scalar avg' endif - call utilities_FFTscalarBackward() + call fftw_mpi_execute_dft_c2r(planScalarBack,scalarField_fourier,scalarField_real) scalarField_real(cells(1)+1:cells1Red*2,:,:) = 0.0_pReal - if (maxval(abs(scalarField_real_ - scalarField_real))>5.0e-15_pReal) error stop 'scalarField' + if (maxval(abs(scalarField_real_ - scalarField_real*wgt))>5.0e-15_pReal) error stop 'scalarField' end subroutine selfTest