avoid use of global variables

src/grid/grid_damage_spectral.f90

@@ -318,16 +318,14 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc)
   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)) &
+    r(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)) &
+             + homogenization_mu_phi(ce)*(phi_lastInc(i,j,k) - phi_current(i,j,k)) &
-                            + mu_ref*phi_current(i,j,k)
+             + mu_ref*phi_current(i,j,k)
   end do; end do; end do
 
-  call utilities_fourierGreenConvolution(K_ref, mu_ref, params%Delta_t)
-
 !--------------------------------------------------------------------------------------------------
 ! constructing residual
-  r = max(min(scalarField_real(1:cells(1),1:cells(2),1:cells3),phi_lastInc),num%residualStiffness) &
+  r = max(min(utilities_GreenConvolution(r, K_ref, mu_ref, params%Delta_t),phi_lastInc),num%residualStiffness) &
     - phi_current
   err_PETSc = 0
 

src/grid/grid_mech_spectral_basic.f90

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

src/grid/grid_mech_spectral_polarisation.f90

@@ -621,7 +621,7 @@ subroutine formResidual(in, FandF_tau, &
 !--------------------------------------------------------------------------------------------------
 ! constructing residual
   r_F_tau = num%beta*F &
-          - utilities_fourierGammaConvolution(r_F_tau,params%rotation_BC%rotate(num%beta*F_aim,active=.true.))
+          - utilities_GammaConvolution(r_F_tau,params%rotation_BC%rotate(num%beta*F_aim,active=.true.))
 
 !--------------------------------------------------------------------------------------------------
 ! evaluate constitutive response

src/grid/grid_thermal_spectral.f90

@@ -341,17 +341,15 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc)
   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)) &
+    r(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)) &
+             + homogenization_mu_T(ce) * (T_lastInc(i,j,k) - T_current(i,j,k)) &
-                            + mu_ref*T_current(i,j,k)
+             + mu_ref*T_current(i,j,k)
   end do; end do; end do
 
-  call utilities_fourierGreenConvolution(K_ref, mu_ref, params%Delta_t)
-
 !--------------------------------------------------------------------------------------------------
 ! constructing residual
   r = T_current &
-    - scalarField_real(1:cells(1),1:cells(2),1:cells3)
+    - utilities_GreenConvolution(r, K_ref, mu_ref, params%Delta_t)
   err_PETSc = 0
 
 end subroutine formResidual

src/grid/spectral_utilities.f90

@@ -116,8 +116,8 @@ module spectral_utilities
   public :: &
     spectral_utilities_init, &
     utilities_updateGamma, &
-    utilities_fourierGammaConvolution, &
+    utilities_GammaConvolution, &
-    utilities_fourierGreenConvolution, &
+    utilities_GreenConvolution, &
     utilities_divergenceRMS, &
     utilities_curlRMS, &
     utilities_fourierScalarGradient, &
@@ -507,7 +507,7 @@ end subroutine utilities_FFTvectorBackward
 !--------------------------------------------------------------------------------------------------
 !> @brief doing convolution gamma_hat * field_real, ensuring that average value = fieldAim
 !--------------------------------------------------------------------------------------------------
-function utilities_fourierGammaConvolution(field, fieldAim) result(gammaField)
+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
@@ -528,8 +528,6 @@ function utilities_fourierGammaConvolution(field, fieldAim) result(gammaField)
   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)
     do j = 1, cells2; do k = 1, cells(3); do i = 1, cells1Red
@@ -592,22 +590,27 @@ function utilities_fourierGammaConvolution(field, fieldAim) result(gammaField)
   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_fourierGammaConvolution
+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
-  call utilities_FFTscalarForward()
+  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) &
@@ -618,7 +621,9 @@ subroutine utilities_fourierGreenConvolution(D_ref, mu_ref, Delta_t)
   !$OMP END PARALLEL DO
   call utilities_FFTscalarBackward()
 
-end subroutine utilities_fourierGreenConvolution
+  greenField = scalarField_real(1:cells(1),1:cells(2),1:cells3)
+
+end function utilities_GreenConvolution
 
 
 !--------------------------------------------------------------------------------------------------