avoid global variables

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
+  r = utilities_fourierGammaConvolution(r,params%rotation_BC%rotate(deltaF_aim,active=.true.))      ! convolution of Gamma and tensorField_fourier
-  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
 
 end subroutine formResidual
 

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

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(3,3,cells(1),cells(2),cells3) :: gammaField
-  real(pReal),             dimension(6,6) :: A, A_inv
 
+  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
 
 
 !--------------------------------------------------------------------------------------------------