internal (private) functions at the end of the module

src/grid_damage_spectral.f90

@@ -170,7 +170,7 @@ function grid_damage_spectral_solution(timeinc,timeinc_old,loadCaseTime) result(
     loadCaseTime                                                                                    !< remaining time of current load case
   integer :: i, j, k, cell
   type(tSolutionState) :: solution
-  PetscInt  ::position
+  PetscInt  :: devNull
   PetscReal :: minDamage, maxDamage, stagNorm, solnNorm
   
   PetscErrorCode :: ierr
@@ -208,8 +208,8 @@ function grid_damage_spectral_solution(timeinc,timeinc_old,loadCaseTime) result(
     call damage_nonlocal_putNonLocalDamage(damage_current(i,j,k),1,cell)
   enddo; enddo; enddo
  
-  call VecMin(solution_vec,position,minDamage,ierr); CHKERRQ(ierr)
+  call VecMin(solution_vec,devNull,minDamage,ierr); CHKERRQ(ierr)
-  call VecMax(solution_vec,position,maxDamage,ierr); CHKERRQ(ierr)
+  call VecMax(solution_vec,devNull,maxDamage,ierr); CHKERRQ(ierr)
   if (solution%converged) &
     write(6,'(/,a)') ' ... nonlocal damage converged .....................................'
   write(6,'(/,a,f8.6,2x,f8.6,2x,f8.6,/)',advance='no') ' Minimum|Maximum|Delta Damage      = ',&

src/grid_mech_spectral_basic.f90

@@ -35,7 +35,9 @@ module grid_mech_spectral_basic
 
 !--------------------------------------------------------------------------------------------------
 ! common pointwise data
-  real(pReal), private, dimension(:,:,:,:,:), allocatable ::  F_lastInc, Fdot
+  real(pReal), private, dimension(:,:,:,:,:), allocatable ::  &
+    F_lastInc, &
+    Fdot
 
 !--------------------------------------------------------------------------------------------------
 ! stress, stiffness and compliance average etc.
@@ -44,9 +46,8 @@ module grid_mech_spectral_basic
     F_aim = math_I3, &                                                                              !< current prescribed deformation gradient
     F_aim_lastInc = math_I3, &                                                                      !< previous average deformation gradient
     P_av = 0.0_pReal                                                                                !< average 1st Piola--Kirchhoff stress
- 
+
-  character(len=1024), private :: incInfo                                                           !< time and increment information
+  character(len=1024), private :: incInfo                                                           !< time and increment information 
- 
   real(pReal), private, dimension(3,3,3,3) :: &
     C_volAvg = 0.0_pReal, &                                                                         !< current volume average stiffness 
     C_volAvgLastInc = 0.0_pReal, &                                                                  !< previous volume average stiffness
@@ -65,6 +66,8 @@ module grid_mech_spectral_basic
     grid_mech_spectral_basic_init, &
     grid_mech_spectral_basic_solution, &
     grid_mech_spectral_basic_forward
+  private :: &
+    formResidual
 
 contains
 
@@ -107,7 +110,8 @@ subroutine grid_mech_spectral_basic_init
     temp33_Real = 0.0_pReal
  
   PetscErrorCode :: ierr
-  PetscScalar, pointer, dimension(:,:,:,:)   ::  F
+  PetscScalar, pointer, dimension(:,:,:,:) :: &
+   F                                                                                                ! pointer to solution data
   PetscInt, dimension(worldsize) :: localK  
   integer :: fileUnit
   character(len=1024) :: rankStr
@@ -152,7 +156,7 @@ subroutine grid_mech_spectral_basic_init
   call DMsetFromOptions(da,ierr); CHKERRQ(ierr)
   call DMsetUp(da,ierr); CHKERRQ(ierr)
   call DMcreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr)                                    ! global solution vector (grid x 9, i.e. every def grad tensor)
-  call DMDASNESsetFunctionLocal(da,INSERT_VALUES,grid_mech_spectral_basic_formResidual,PETSC_NULL_SNES,ierr)! residual vector of same shape as solution vector
+  call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,ierr)                 ! residual vector of same shape as solution vector
   CHKERRQ(ierr) 
   call SNESsetConvergenceTest(snes,grid_mech_spectral_basic_converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr)! specify custom convergence check function "_converged"
   CHKERRQ(ierr)
@@ -160,7 +164,7 @@ subroutine grid_mech_spectral_basic_init
 
 !--------------------------------------------------------------------------------------------------
 ! init fields                 
-  call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)                                    ! get the data out of PETSc to work with
+  call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)                                   ! places pointer on PETSc data
  
   restart: if (restartInc > 0) then                                                     
     if (iand(debug_level(debug_spectral),debug_spectralRestart) /= 0) then
@@ -196,11 +200,10 @@ subroutine grid_mech_spectral_basic_init
                                       reshape(F,shape(F_lastInc)), &                                ! target F
                                       0.0_pReal, &                                                  ! time increment
                                       math_I3)                                                      ! no rotation of boundary condition
-  call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)                                ! write data back to PETSc
+  call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)                                ! deassociate pointer
-                                                                                                    ! QUESTION: why not writing back right after reading (l.189)?
  
   restartRead: if (restartInc > 0) then
-    if (iand(debug_level(debug_spectral),debug_spectralRestart) /= 0) &
+    if (iand(debug_level(debug_spectral),debug_spectralRestart) /= 0 .and. worldrank == 0) &
       write(6,'(/,a,'//IO_intOut(restartInc)//',a)') &
       'reading more values of increment ', restartInc, ' from file'
     flush(6)
@@ -218,7 +221,7 @@ end subroutine grid_mech_spectral_basic_init
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief solution for the Basic scheme with internal iterations
+!> @brief solution for the basic scheme with internal iterations
 !--------------------------------------------------------------------------------------------------
 function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution)
   use numerics, only: &
@@ -245,7 +248,6 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_
   real(pReal), dimension(3,3), intent(in) :: rotation_BC
   type(tSolutionState)                    :: &
     solution
- 
 !--------------------------------------------------------------------------------------------------
 ! PETSc Data
   PetscErrorCode :: ierr
@@ -283,97 +285,6 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_
 end function grid_mech_spectral_basic_solution
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief forms the basic residual vector
-!--------------------------------------------------------------------------------------------------
-subroutine grid_mech_spectral_basic_formResidual(in, F, &
-                              residuum, dummy, ierr)
-  use numerics, only: &
-    itmax, &
-    itmin
-  use mesh, only: &
-    grid, &
-    grid3
-  use math, only: &
-    math_rotate_backward33, &
-    math_mul3333xx33
-  use debug, only: &
-    debug_level, &
-    debug_spectral, &
-    debug_spectralRotation
-  use spectral_utilities, only: &
-    tensorField_real, &
-    utilities_FFTtensorForward, &
-    utilities_fourierGammaConvolution, &
-    utilities_FFTtensorBackward, &
-    utilities_constitutiveResponse, &
-    utilities_divergenceRMS
-  use IO, only: &
-    IO_intOut 
-  use FEsolving, only: &
-    terminallyIll
-
-  implicit none
-  DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in                                              !< DMDA info (needs to be named "in" for macros like XRANGE to work)
-  PetscScalar, dimension(3,3,XG_RANGE,YG_RANGE,ZG_RANGE), &
-    intent(in)  :: F                                                                                !< deformation gradient field
-  PetscScalar, dimension(3,3,X_RANGE,Y_RANGE,Z_RANGE), &
-    intent(out) :: residuum                                                                         !< residuum field
-  real(pReal),  dimension(3,3) :: &
-    deltaF_aim
-  PetscInt :: &
-    PETScIter, &
-    nfuncs
-  PetscObject :: dummy
-  PetscErrorCode :: ierr
-
-  call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr)
-  call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr)
-
-  if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1                                              ! new increment
-!--------------------------------------------------------------------------------------------------
-! begin of new iteration
-  newIteration: if (totalIter <= PETScIter) then
-    totalIter = totalIter + 1
-    write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') &
-            trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
-    if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
-      write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
-              ' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC))
-    write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
-              ' deformation gradient aim       =', transpose(F_aim)
-    flush(6)
-  endif newIteration
-
-!--------------------------------------------------------------------------------------------------
-! evaluate constitutive response
-  call Utilities_constitutiveResponse(residuum, &                                                   ! "residuum" gets field of first PK stress (to save memory)
-                                      P_av,C_volAvg,C_minMaxAvg, &
-                                      F,params%timeinc,params%rotation_BC)
-  call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)
-  
-!--------------------------------------------------------------------------------------------------
-! stress BC handling
-  deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC)
-  F_aim = F_aim - deltaF_aim
-  err_BC = maxval(abs(params%stress_mask * (P_av - params%stress_BC)))                              ! mask = 0.0 when no stress bc
-
-!--------------------------------------------------------------------------------------------------
-! updated deformation gradient using fix point algorithm of basic scheme
-  tensorField_real = 0.0_pReal
-  tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residuum                                  ! 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(math_rotate_backward33(deltaF_aim,params%rotation_BC))     ! convolution of Gamma and tensorField_fourier, with arg 
-  call utilities_FFTtensorBackward                                                                  ! FFT backward of global tensorField_fourier
- 
-!--------------------------------------------------------------------------------------------------
-! constructing residual
-  residuum = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3)                                   ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too
-
-end subroutine grid_mech_spectral_basic_formResidual
-
-
 !--------------------------------------------------------------------------------------------------
 !> @brief convergence check
 !--------------------------------------------------------------------------------------------------
@@ -542,4 +453,96 @@ subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTi
   
 end subroutine grid_mech_spectral_basic_forward
 
+
+!--------------------------------------------------------------------------------------------------
+!> @brief forms the basic residual vector
+!--------------------------------------------------------------------------------------------------
+subroutine formResidual(in, F, &
+                        residuum, dummy, ierr)
+  use numerics, only: &
+    itmax, &
+    itmin
+  use mesh, only: &
+    grid, &
+    grid3
+  use math, only: &
+    math_rotate_backward33, &
+    math_mul3333xx33
+  use debug, only: &
+    debug_level, &
+    debug_spectral, &
+    debug_spectralRotation
+  use spectral_utilities, only: &
+    tensorField_real, &
+    utilities_FFTtensorForward, &
+    utilities_fourierGammaConvolution, &
+    utilities_FFTtensorBackward, &
+    utilities_constitutiveResponse, &
+    utilities_divergenceRMS
+  use IO, only: &
+    IO_intOut 
+  use FEsolving, only: &
+    terminallyIll
+
+  implicit none
+  DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in                                              !< DMDA info (needs to be named "in" for macros like XRANGE to work)
+  PetscScalar, dimension(3,3,XG_RANGE,YG_RANGE,ZG_RANGE), &
+    intent(in) :: F                                                                                 !< deformation gradient field
+  PetscScalar, dimension(3,3,X_RANGE,Y_RANGE,Z_RANGE), &
+    intent(out) :: residuum                                                                         !< residuum field
+  real(pReal),  dimension(3,3) :: &
+    deltaF_aim
+  PetscInt :: &
+    PETScIter, &
+    nfuncs
+  PetscObject :: dummy
+  PetscErrorCode :: ierr
+
+  call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr)
+  call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr)
+
+  if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1                                              ! new increment
+!--------------------------------------------------------------------------------------------------
+! begin of new iteration
+  newIteration: if (totalIter <= PETScIter) then
+    totalIter = totalIter + 1
+    write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') &
+            trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
+    if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
+      write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
+              ' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC))
+    write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
+              ' deformation gradient aim       =', transpose(F_aim)
+    flush(6)
+  endif newIteration
+
+!--------------------------------------------------------------------------------------------------
+! evaluate constitutive response
+  call Utilities_constitutiveResponse(residuum, &                                                   ! "residuum" gets field of first PK stress (to save memory)
+                                      P_av,C_volAvg,C_minMaxAvg, &
+                                      F,params%timeinc,params%rotation_BC)
+  call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)
+  
+!--------------------------------------------------------------------------------------------------
+! stress BC handling
+  deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC)
+  F_aim = F_aim - deltaF_aim
+  err_BC = maxval(abs(params%stress_mask * (P_av - params%stress_BC)))                              ! mask = 0.0 when no stress bc
+
+!--------------------------------------------------------------------------------------------------
+! updated deformation gradient using fix point algorithm of basic scheme
+  tensorField_real = 0.0_pReal
+  tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residuum                                  ! 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(math_rotate_backward33(deltaF_aim,params%rotation_BC))     ! convolution of Gamma and tensorField_fourier, with arg 
+  call utilities_FFTtensorBackward                                                                  ! FFT backward of global tensorField_fourier
+ 
+!--------------------------------------------------------------------------------------------------
+! constructing residual
+  residuum = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3)                                   ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too
+
+end subroutine formResidual
+
+
 end module grid_mech_spectral_basic

src/grid_mech_spectral_polarisation.f90

@@ -2,7 +2,7 @@
 !> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
 !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
 !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
-!> @brief Polarisation scheme solver
+!> @brief Grid solver for mechanics: Spectral Polarisation
 !--------------------------------------------------------------------------------------------------
 module grid_mech_spectral_polarisation
 #include <petsc/finclude/petscsnes.h>
@@ -10,7 +10,6 @@ module grid_mech_spectral_polarisation
  use PETScdmda
  use PETScsnes
  use prec, only: & 
-   pInt, &
    pReal
  use math, only: &
    math_I3
@@ -51,7 +50,8 @@ module grid_mech_spectral_polarisation
    F_av = 0.0_pReal, &                                                                              !< average incompatible def grad field
    P_av = 0.0_pReal, &                                                                              !< average 1st Piola--Kirchhoff stress
    P_avLastEval = 0.0_pReal                                                                         !< average 1st Piola--Kirchhoff stress last call of CPFEM_general
-   character(len=1024), private :: incInfo                                                          !< time and increment information
+ 
+ character(len=1024), private :: incInfo                                                          !< time and increment information
  real(pReal), private, dimension(3,3,3,3) :: &
    C_volAvg = 0.0_pReal, &                                                                          !< current volume average stiffness 
    C_volAvgLastInc = 0.0_pReal, &                                                                   !< previous volume average stiffness
@@ -66,13 +66,15 @@ module grid_mech_spectral_polarisation
    err_curl, &                                                                                      !< RMS of curl of F
    err_div                                                                                          !< RMS of div of P
   
- integer(pInt), private :: &
+ integer, private :: &
-   totalIter = 0_pInt                                                                               !< total iteration in current increment
+   totalIter = 0                                                                                    !< total iteration in current increment
  
  public :: &
    grid_mech_spectral_polarisation_init, &
    grid_mech_spectral_polarisation_solution, &
    grid_mech_spectral_polarisation_forward
+ private :: &
+   formResidual
 
 contains
 
@@ -99,9 +101,9 @@ subroutine grid_mech_spectral_polarisation_init
  use DAMASK_interface, only: &
    getSolverJobName
  use spectral_utilities, only: &
-   Utilities_constitutiveResponse, &
+    utilities_constitutiveResponse, &
-   Utilities_updateGamma, &
+    utilities_updateGamma, &
-   Utilities_updateIPcoords, &
+    utilities_updateIPcoords, &
    wgt
  use mesh, only: &
    grid, &
@@ -123,7 +125,7 @@ subroutine grid_mech_spectral_polarisation_init
  integer :: fileUnit
  character(len=1024) :: rankStr
  
- write(6,'(/,a)') ' <<<+-  DAMASK_spectral_solverPolarisation init  -+>>>'
+ write(6,'(/,a)') ' <<<+-  grid_mech_spectral_polarisation init  -+>>>'
 
  write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
  write(6,'(a)')   ' https://doi.org/10.1016/j.ijplas.2014.02.006'
@@ -189,7 +191,6 @@ subroutine grid_mech_spectral_polarisation_init
    read(fileUnit) F; close (fileUnit)
    fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr))
    read(fileUnit) F_lastInc; close (fileUnit)
-
    fileUnit = IO_open_jobFile_binary('F_tau'//trim(rankStr))
    read(fileUnit) F_tau; close (fileUnit)
    fileUnit = IO_open_jobFile_binary('F_tau_lastInc'//trim(rankStr))
@@ -197,11 +198,11 @@ subroutine grid_mech_spectral_polarisation_init
 
    F_aim         = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3])                         ! average of F
    call MPI_Allreduce(MPI_IN_PLACE,F_aim,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
-   if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='F_aim')
+   if(ierr /=0) call IO_error(894, ext_msg='F_aim')
    F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt                                 ! average of F_lastInc 
    call MPI_Allreduce(MPI_IN_PLACE,F_aim_lastInc,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
-   if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='F_aim_lastInc')
+   if(ierr /=0) call IO_error(894, ext_msg='F_aim_lastInc')
- elseif (restartInc == 0_pInt) then restart
+ elseif (restartInc == 0) then restart
    F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3)                          ! initialize to identity
    F = reshape(F_lastInc,[9,grid(1),grid(2),grid3])
    F_tau = 2.0_pReal*F
@@ -214,12 +215,10 @@ subroutine grid_mech_spectral_polarisation_init
                                      reshape(F,shape(F_lastInc)), &                                 ! target F
                                      0.0_pReal, &                                                   ! time increment
                                      math_I3)                                                       ! no rotation of boundary condition
- nullify(F)
+ call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)                         ! deassociate pointer
- nullify(F_tau)
- call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)                         ! write data back to PETSc
 
- restartRead: if (restartInc > 0_pInt) then
+ restartRead: if (restartInc > 0) then
-   if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) &
+   if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0) &
      write(6,'(/,a,'//IO_intOut(restartInc)//',a)') &
      'reading more values of increment ', restartInc, ' from file'
    flush(6)
@@ -250,8 +249,8 @@ function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old,
    math_invSym3333
  use spectral_utilities, only: &
    tBoundaryCondition, &
-   Utilities_maskedCompliance, &
+    utilities_maskedCompliance, &
-   Utilities_updateGamma
+    utilities_updateGamma
  use FEsolving, only: &
    restartWrite, &
    terminallyIll
@@ -287,7 +286,7 @@ function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old,
  endif  
 
 !--------------------------------------------------------------------------------------------------
-! set module wide availabe data
+! set module wide available data 
  params%stress_mask = stress_BC%maskFloat
  params%stress_BC   = stress_BC%values
  params%rotation_BC = rotation_BC
@@ -306,161 +305,10 @@ function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old,
  solution%iterationsNeeded = totalIter
  solution%termIll = terminallyIll
  terminallyIll = .false.
- if (reason == -4) call IO_error(893_pInt)                                                         ! MPI error
 
 end function grid_mech_spectral_polarisation_solution
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief forms the Polarisation residual vector
-!--------------------------------------------------------------------------------------------------
-subroutine formResidual(in, &                                                         ! DMDA info (needs to be named "in" for XRANGE, etc. macros to work)
-                                     FandF_tau, &                                                  ! defgrad fields on grid
-                                     residuum, &                                                   ! residuum fields on grid
-                                     dummy, &
-                                     ierr)
- use numerics, only: &
-   itmax, &
-   itmin, &
-   polarAlpha, &
-   polarBeta
- use mesh, only: &
-   grid, &
-   grid3
- use IO, only: &
-   IO_intOut
- use math, only: &
-   math_rotate_backward33, &
-   math_mul3333xx33, &
-   math_invSym3333, &
-   math_mul33x33
- use debug, only: &
-   debug_level, &
-   debug_spectral, &
-   debug_spectralRotation
- use spectral_utilities, only: &
-   wgt, &
-   tensorField_real, &
-   utilities_FFTtensorForward, &
-   utilities_fourierGammaConvolution, &
-   utilities_FFTtensorBackward, &
-   Utilities_constitutiveResponse, &
-   Utilities_divergenceRMS, &
-   Utilities_curlRMS
- use homogenization, only: &
-   materialpoint_dPdF
- use FEsolving, only: &
-   terminallyIll
-
- implicit none
- DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in
- PetscScalar, &
-   target, dimension(3,3,2, XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: FandF_tau
- PetscScalar, &
-   target, dimension(3,3,2,  X_RANGE, Y_RANGE, Z_RANGE), intent(out) :: residuum
- PetscScalar, pointer, dimension(:,:,:,:,:) :: &
-   F, &
-   F_tau, &
-   residual_F, &
-   residual_F_tau
- PetscInt :: &
-   PETScIter, &
-   nfuncs
- PetscObject :: dummy
- PetscErrorCode :: ierr
- integer(pInt) :: &
-   i, j, k, e
-
- F                 => FandF_tau(1:3,1:3,1,&
-                                XG_RANGE,YG_RANGE,ZG_RANGE)
- F_tau             => FandF_tau(1:3,1:3,2,&
-                                XG_RANGE,YG_RANGE,ZG_RANGE)
- residual_F        => residuum(1:3,1:3,1,&
-                                X_RANGE, Y_RANGE, Z_RANGE)
- residual_F_tau    => residuum(1:3,1:3,2,&
-                                X_RANGE, Y_RANGE, Z_RANGE)
-
- F_av = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt
- call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
- 
- call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr)
- call SNESGetIterationNumber(snes,PETScIter,ierr);  CHKERRQ(ierr)
-
- if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1_pInt                                            ! new increment
-!--------------------------------------------------------------------------------------------------
-! begin of new iteration
- newIteration: if (totalIter <= PETScIter) then
-   totalIter = totalIter + 1_pInt
-   write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') &
-           trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
-   if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
-     write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
-             ' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC))
-   write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
-             ' deformation gradient aim       =', transpose(F_aim)
-   flush(6)
- endif newIteration
-
-!--------------------------------------------------------------------------------------------------
-! 
- tensorField_real = 0.0_pReal
- do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid(1)
-   tensorField_real(1:3,1:3,i,j,k) = &
-     polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -&
-     polarAlpha*math_mul33x33(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))
- enddo; enddo; enddo
- 
-!--------------------------------------------------------------------------------------------------
-! doing convolution in Fourier space 
- call utilities_FFTtensorForward()
- call utilities_fourierGammaConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC)) 
- call utilities_FFTtensorBackward()
-
-!--------------------------------------------------------------------------------------------------
-! constructing residual                         
- residual_F_tau = polarBeta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3)
-
-!--------------------------------------------------------------------------------------------------
-! evaluate constitutive response
- P_avLastEval = P_av
- call Utilities_constitutiveResponse(residual_F,P_av,C_volAvg,C_minMaxAvg, &
-                                     F - residual_F_tau/polarBeta,params%timeinc,params%rotation_BC)
- call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) 
-
-!--------------------------------------------------------------------------------------------------
-! calculate divergence
- tensorField_real = 0.0_pReal
- tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F                                   !< stress field in disguise
- call utilities_FFTtensorForward()
- err_div = Utilities_divergenceRMS()                                                                  !< root mean squared error in divergence of stress
- 
-!--------------------------------------------------------------------------------------------------
-! constructing residual
- e = 0_pInt
- do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid(1)
-   e = e + 1_pInt
-   residual_F(1:3,1:3,i,j,k) = &
-     math_mul3333xx33(math_invSym3333(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,e) + C_scale), &
-                      residual_F(1:3,1:3,i,j,k) - math_mul33x33(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))) &
-                      + residual_F_tau(1:3,1:3,i,j,k)
- enddo; enddo; enddo
- 
-!--------------------------------------------------------------------------------------------------
-! calculating curl
- tensorField_real = 0.0_pReal
- tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F
- call utilities_FFTtensorForward()
- err_curl = Utilities_curlRMS()
-
- nullify(F)
- nullify(F_tau)
- nullify(residual_F)
- nullify(residual_F_tau)
-end subroutine formResidual
-
-
 !--------------------------------------------------------------------------------------------------
 !> @brief convergence check
 !--------------------------------------------------------------------------------------------------
@@ -483,9 +331,9 @@ subroutine grid_mech_spectral_polarisation_converged(snes_local,PETScIter,xnorm,
  SNES :: snes_local
  PetscInt :: PETScIter
  PetscReal :: &
-   xnorm, &
+   xnorm, &                                                                                        ! not used
-   snorm, &
+   snorm, &                                                                                        ! not used
-   fnorm
+   fnorm                                                                                           ! not used
  SNESConvergedReason :: reason
  PetscObject :: dummy
  PetscErrorCode :: ierr
@@ -552,9 +400,9 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa
   use CPFEM2, only: &
     CPFEM_age
   use spectral_utilities, only: &
-    Utilities_calculateRate, &
+    utilities_calculateRate, &
-    Utilities_forwardField, &
+    utilities_forwardField, &
-    Utilities_updateIPcoords, &
+    utilities_updateIPcoords, &
     tBoundaryCondition, &
     cutBack
   use IO, only: &
@@ -576,7 +424,7 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa
     rotation_BC
   PetscErrorCode :: ierr
   PetscScalar, dimension(:,:,:,:), pointer :: FandF_tau, F, F_tau
-  integer(pInt) :: i, j, k
+  integer :: i, j, k
   real(pReal), dimension(3,3) :: F_lambda33
 
   integer :: fileUnit
@@ -664,7 +512,7 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa
     F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot), &
                     [9,grid(1),grid(2),grid3])                                                        ! does not have any average value as boundary condition
   else
-   do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid(1)
+   do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1)
       F_lambda33 = reshape(F_tau(1:9,i,j,k)-F(1:9,i,j,k),[3,3])
       F_lambda33 = math_mul3333xx33(S_scale,math_mul33x33(F_lambda33, &
                                   math_mul3333xx33(C_scale,&
@@ -681,4 +529,148 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa
 
 end subroutine grid_mech_spectral_polarisation_forward
 
+
+!--------------------------------------------------------------------------------------------------
+!> @brief forms the polarisation residual vector
+!--------------------------------------------------------------------------------------------------
+subroutine formResidual(in, FandF_tau, &
+                        residuum, dummy,ierr)
+ use numerics, only: &
+   itmax, &
+   itmin, &
+   polarAlpha, &
+   polarBeta
+ use mesh, only: &
+   grid, &
+   grid3
+ use IO, only: &
+   IO_intOut
+ use math, only: &
+   math_rotate_backward33, &
+   math_mul3333xx33, &
+   math_invSym3333, &
+   math_mul33x33
+ use debug, only: &
+   debug_level, &
+   debug_spectral, &
+   debug_spectralRotation
+ use spectral_utilities, only: &
+   wgt, &
+   tensorField_real, &
+   utilities_FFTtensorForward, &
+   utilities_fourierGammaConvolution, &
+   utilities_FFTtensorBackward, &
+   Utilities_constitutiveResponse, &
+   Utilities_divergenceRMS, &
+   Utilities_curlRMS
+ use homogenization, only: &
+   materialpoint_dPdF
+ use FEsolving, only: &
+   terminallyIll
+
+ implicit none
+ DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in                                               !< DMDA info (needs to be named "in" for macros like XRANGE to work)
+ PetscScalar, dimension(3,3,2, XG_RANGE,YG_RANGE,ZG_RANGE), &
+   target, intent(in) :: FandF_tau
+ PetscScalar, dimension(3,3,2,  X_RANGE, Y_RANGE, Z_RANGE),&
+   target,  intent(out) :: residuum                                                                 !< residuum field
+ PetscScalar, pointer, dimension(:,:,:,:,:) :: &
+   F, &
+   F_tau, &
+   residual_F, &
+   residual_F_tau
+ PetscInt :: &
+   PETScIter, &
+   nfuncs
+ PetscObject :: dummy
+ PetscErrorCode :: ierr
+ integer :: &
+   i, j, k, e
+
+ F                 => FandF_tau(1:3,1:3,1,&
+                                XG_RANGE,YG_RANGE,ZG_RANGE)
+ F_tau             => FandF_tau(1:3,1:3,2,&
+                                XG_RANGE,YG_RANGE,ZG_RANGE)
+ residual_F        => residuum(1:3,1:3,1,&
+                                X_RANGE, Y_RANGE, Z_RANGE)
+ residual_F_tau    => residuum(1:3,1:3,2,&
+                                X_RANGE, Y_RANGE, Z_RANGE)
+
+ F_av = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt
+ call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
+ 
+ call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr)
+ call SNESGetIterationNumber(snes,PETScIter,ierr);  CHKERRQ(ierr)
+
+ if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1                                                ! new increment
+!--------------------------------------------------------------------------------------------------
+! begin of new iteration
+ newIteration: if (totalIter <= PETScIter) then
+   totalIter = totalIter + 1
+   write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') &
+           trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
+   if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
+     write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
+             ' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC))
+   write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
+             ' deformation gradient aim       =', transpose(F_aim)
+   flush(6)
+ endif newIteration
+
+!--------------------------------------------------------------------------------------------------
+! 
+ tensorField_real = 0.0_pReal
+ do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1)
+   tensorField_real(1:3,1:3,i,j,k) = &
+     polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -&
+     polarAlpha*math_mul33x33(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))
+ enddo; enddo; enddo
+ 
+!--------------------------------------------------------------------------------------------------
+! doing convolution in Fourier space 
+ call utilities_FFTtensorForward()
+ call utilities_fourierGammaConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC)) 
+ call utilities_FFTtensorBackward()
+
+!--------------------------------------------------------------------------------------------------
+! constructing residual                         
+ residual_F_tau = polarBeta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3)
+
+!--------------------------------------------------------------------------------------------------
+! evaluate constitutive response
+ P_avLastEval = P_av
+ call Utilities_constitutiveResponse(residual_F, &                                                    ! "residuum" gets field of first PK stress (to save memory)
+                                     P_av,C_volAvg,C_minMaxAvg, &
+                                     F - residual_F_tau/polarBeta,params%timeinc,params%rotation_BC)
+ call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) 
+
+!--------------------------------------------------------------------------------------------------
+! calculate divergence
+ tensorField_real = 0.0_pReal
+ tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_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, grid3; do j = 1, grid(2); do i = 1, grid(1)
+   e = e + 1
+   residual_F(1:3,1:3,i,j,k) = &
+     math_mul3333xx33(math_invSym3333(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,e) + C_scale), &
+                      residual_F(1:3,1:3,i,j,k) - math_mul33x33(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))) &
+                      + residual_F_tau(1:3,1:3,i,j,k)
+ enddo; enddo; enddo
+ 
+!--------------------------------------------------------------------------------------------------
+! calculating curl
+ tensorField_real = 0.0_pReal
+ tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F
+ call utilities_FFTtensorForward()
+ err_curl = Utilities_curlRMS()
+
+end subroutine formResidual
+
 end module grid_mech_spectral_polarisation