diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index a065633de..1ec6d334c 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) @@ -544,12 +546,7 @@ subroutine formResidual(in, F, & F_aim = F_aim - deltaF_aim err_BC = maxval(abs(merge(.0_pReal,P_av - P_aim,params%stress_mask))) - 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_fourierGammaConvolution(params%rotation_BC%rotate(deltaF_aim,active=.true.)) ! convolution of Gamma and 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_fourierGammaConvolution(r,params%rotation_BC%rotate(deltaF_aim,active=.true.)) ! convolution of Gamma and tensorField_fourier end subroutine formResidual diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index fea34d7d0..235bb873e 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) @@ -612,17 +612,16 @@ subroutine formResidual(in, FandF_tau, & !-------------------------------------------------------------------------------------------------- ! 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 - call utilities_fourierGammaConvolution(params%rotation_BC%rotate(num%beta*F_aim,active=.true.)) - !-------------------------------------------------------------------------------------------------- ! constructing residual - r_F_tau = num%beta*F - tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) + r_F_tau = num%beta*F & + - utilities_fourierGammaConvolution(r_F_tau,params%rotation_BC%rotate(num%beta*F_aim,active=.true.)) !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 858fc20db..9e5282de0 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -42,7 +42,7 @@ 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), 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 @@ -505,12 +505,14 @@ end subroutine utilities_FFTvectorBackward !-------------------------------------------------------------------------------------------------- !> @brief doing convolution gamma_hat * field_real, ensuring that average value = fieldAim !-------------------------------------------------------------------------------------------------- -subroutine utilities_fourierGammaConvolution(fieldAim) +function utilities_fourierGammaConvolution(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 @@ -520,8 +522,11 @@ subroutine utilities_fourierGammaConvolution(fieldAim) print'(/,1x,a)', '... doing gamma convolution ...............................................' flush(IO_STDOUT) - call utilities_FFTtensorForward() -!-------------------------------------------------------------------------------------------------- + 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) + + !-------------------------------------------------------------------------------------------------- ! do the actual spectral method calculation (mechanical equilibrium) memoryEfficient: if (num%memory_efficient) then !$OMP PARALLEL DO PRIVATE(l,m,n,o,temp33_cmplx,xiDyad_cmplx,A,A_inv,err,gamma_hat) @@ -582,9 +587,10 @@ subroutine utilities_fourierGammaConvolution(fieldAim) end if memoryEfficient if (cells3Offset == 0) tensorField_fourier(1:3,1:3,1,1,1) = cmplx(fieldAim/wgt,0.0_pReal,pReal) - call utilities_FFTtensorBackward() + 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)*wgt -end subroutine utilities_fourierGammaConvolution +end function utilities_fourierGammaConvolution !--------------------------------------------------------------------------------------------------