separated forwarding and solution subroutines for better control at the load step looping level

code/DAMASK_spectral_driver.f90

@@ -433,6 +433,37 @@ program DAMASK_spectral_Driver
          stepFraction = stepFraction + 1_pInt 
          remainingLoadCaseTime = time0 - time + loadCases(currentLoadCase)%time + timeInc
 !--------------------------------------------------------------------------------------------------
+! forward solution 
+         select case(myspectralsolver)
+           case (DAMASK_spectral_SolverBasic_label)
+#ifdef PETSc
+           case (DAMASK_spectral_SolverBasicPETSC_label)
+             call BasicPETSC_forward (&
+                 guess,timeinc,timeIncOld,remainingLoadCaseTime, &
+                 P_BC               = loadCases(currentLoadCase)%P, &
+                 F_BC               = loadCases(currentLoadCase)%deformation, &
+                 rotation_BC        = loadCases(currentLoadCase)%rotation)
+           
+           case (DAMASK_spectral_SolverAL_label)
+             call AL_forward (&
+                 guess,timeinc,timeIncOld,remainingLoadCaseTime, &
+                 P_BC               = loadCases(currentLoadCase)%P, &
+                 F_BC               = loadCases(currentLoadCase)%deformation, &
+                 rotation_BC        = loadCases(currentLoadCase)%rotation)
+           
+           case (DAMASK_spectral_SolverPolarisation_label)
+             call Polarisation_forward (&
+                 guess,timeinc,timeIncOld,remainingLoadCaseTime, &
+                 P_BC               = loadCases(currentLoadCase)%P, &
+                 F_BC               = loadCases(currentLoadCase)%deformation, &
+                 rotation_BC        = loadCases(currentLoadCase)%rotation)
+#endif
+         end select 
+
+          end select
+         enddo
+           
+!--------------------------------------------------------------------------------------------------
 ! report begin of new increment
          write(6,'(/,a)') ' ###########################################################################'
          write(6,'(1x,a,es12.5'//&

code/DAMASK_spectral_solverAL.f90

@@ -100,6 +100,7 @@ module DAMASK_spectral_solverAL
  public :: &
    AL_init, &
    AL_solution, &
+   AL_forward, &
    AL_destroy
  external :: &
    VecDestroy, &
@@ -277,25 +278,15 @@ type(tSolutionState) function &
  use numerics, only: &
    update_gamma
  use math, only: &
-   math_mul33x33 ,&
+   math_invSym3333
-   math_mul3333xx33, &
-   math_rotate_backward33, &
-   math_invSym3333, &
-   math_transpose33
-use mesh, only: &
-   mesh_ipCoordinates, &
-   mesh_deformedCoordsFFT
  use IO, only: &
    IO_write_jobRealFile
  use DAMASK_spectral_Utilities, only: &
    grid, &
    geomSize, &
    tBoundaryCondition, &
-   Utilities_forwardField, &
-   Utilities_calculateRate, &
    Utilities_maskedCompliance, &
-   Utilities_updateGamma, &
+   Utilities_updateGamma
-   cutBack
  use FEsolving, only: &
    restartWrite, &
    terminallyIll
@@ -327,16 +318,13 @@ use mesh, only: &
  PetscErrorCode :: ierr   
  SNESConvergedReason :: reason
 
- integer(pInt) :: i, j, k
- real(pReal), dimension(3,3) :: F_lambda33
-
  incInfo = incInfoIn
- call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr)
- F => xx_psc(0:8,:,:,:)
- F_lambda => xx_psc(9:17,:,:,:)
  
 !--------------------------------------------------------------------------------------------------
 ! restart information for spectral solver
+ call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr)
+ F => xx_psc(0:8,:,:,:)
+ F_lambda => xx_psc(9:17,:,:,:)
  if (restartWrite) then
    write(6,'(/,a)') ' writing converged results for restart'
    flush(6)
@@ -365,61 +353,10 @@ use mesh, only: &
    write (777,rec=1) C_volAvgLastInc
    close(777)
  endif 
+ call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr)
+ 
  AL_solution%converged =.false.
 
- if (cutBack) then 
-   F_aim = F_aim_lastInc
-   F_lambda= reshape(F_lambda_lastInc,[9,grid(1),grid(2),grid(3)]) 
-   F    = reshape(F_lastInc,    [9,grid(1),grid(2),grid(3)]) 
-   C_volAvg = C_volAvgLastInc
- else
-   C_volAvgLastInc = C_volAvg
-!--------------------------------------------------------------------------------------------------
-! calculate rate for aim
-   if (F_BC%myType=='l') then                                                                       ! calculate f_aimDot from given L and current F
-     f_aimDot = F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim)
-   elseif(F_BC%myType=='fdot') then                                                                 ! f_aimDot is prescribed
-     f_aimDot = F_BC%maskFloat * F_BC%values
-   elseif(F_BC%myType=='f') then                                                                    ! aim at end of load case is prescribed
-     f_aimDot = F_BC%maskFloat * (F_BC%values -F_aim)/loadCaseTime
-  endif
-   if (guess) f_aimDot  = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old
-   F_aim_lastInc = F_aim
-
-!--------------------------------------------------------------------------------------------------
-! update coordinates and rate and forward last inc
-   mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(&
-                                            F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)])
-   Fdot =  Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
-                  timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid(3)]))
-   F_lambdaDot =  Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
-                  timeinc_old,guess,F_lambda_lastInc,reshape(F_lambda,[3,3,grid(1),grid(2),grid(3)]))  
-   F_lastInc2 = F_lastInc                
-   F_lastInc     = reshape(F,       [3,3,grid(1),grid(2),grid(3)])
-   F_lambda_lastInc = reshape(F_lambda,[3,3,grid(1),grid(2),grid(3)])
- endif
- F_aim = F_aim + f_aimDot * timeinc
-
-!--------------------------------------------------------------------------------------------------
-! update local deformation gradient
- F     = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, &                                   ! ensure that it matches rotated F_aim
-                               math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid(3)])
- F_lambda = reshape(Utilities_forwardField(timeinc,F_lambda_lastInc,F_lambdadot), &
-                                                                       [9,grid(1),grid(2),grid(3)]) ! does not have any average value as boundary condition
- if (.not. guess) then                                                                              ! large strain forwarding
-   do k = 1_pInt, grid(3); do j = 1_pInt, grid(2); do i = 1_pInt, grid(1)
-      F_lambda33 = reshape(F_lambda(:,i,j,k),[3,3])
-      F_lambda33 = math_mul3333xx33(S_scale,math_mul33x33(F_lambda33, &
-                                  math_mul3333xx33(C_scale,&
-                                                   math_mul33x33(math_transpose33(F_lambda33),&
-                                                                 F_lambda33) -math_I3))*0.5_pReal)&
-                              + math_I3
-      F_lambda(:,i,j,k) = reshape(F_lambda33,[9])
-   enddo; enddo; enddo
- endif
- call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr)
- CHKERRQ(ierr)
-
 !--------------------------------------------------------------------------------------------------
 ! update stiffness (and gamma operator)
  S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C_volAvg)
@@ -696,6 +633,101 @@ subroutine AL_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr
 
 end subroutine AL_converged
 
+!--------------------------------------------------------------------------------------------------
+!> @brief forwarding routine
+!--------------------------------------------------------------------------------------------------
+subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_BC)
+ use math, only: &
+   math_mul33x33, &
+   math_mul3333xx33, &
+   math_transpose33, &
+   math_rotate_backward33
+ use DAMASK_spectral_Utilities, only: &
+   grid, &
+   geomSize, &
+   Utilities_calculateRate, &
+   Utilities_forwardField, &
+   tBoundaryCondition, &
+   cutBack
+ use mesh, only: &
+   mesh_ipCoordinates,&
+   mesh_deformedCoordsFFT
+
+ implicit none
+#include <finclude/petscdmda.h90>
+ real(pReal), intent(in) :: &
+   timeinc_old, &
+   timeinc, &
+   loadCaseTime                                                                                     !< remaining time of current load case
+ type(tBoundaryCondition),      intent(in) :: &
+   P_BC, &
+   F_BC
+ real(pReal), dimension(3,3), intent(in) :: rotation_BC
+ logical, intent(in) :: &
+   guess
+ PetscErrorCode :: ierr
+ PetscScalar, dimension(:,:,:,:), pointer :: xx_psc, F, F_lambda
+ integer(pInt) :: i, j, k
+ real(pReal), dimension(3,3) :: F_lambda33
+
+!--------------------------------------------------------------------------------------------------
+! update coordinates and rate and forward last inc
+ call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr)
+ F => xx_psc(0:8,:,:,:)
+ F_lambda => xx_psc(9:17,:,:,:) 
+ if (cutBack) then 
+   F_aim = F_aim_lastInc
+   F_lambda= reshape(F_lambda_lastInc,[9,grid(1),grid(2),grid(3)]) 
+   F    = reshape(F_lastInc,    [9,grid(1),grid(2),grid(3)]) 
+   C_volAvg = C_volAvgLastInc
+ else
+!--------------------------------------------------------------------------------------------------
+! calculate rate for aim
+   C_volAvgLastInc = C_volAvg
+   if (F_BC%myType=='l') then                                                                       ! calculate f_aimDot from given L and current F
+     f_aimDot = F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim)
+   elseif(F_BC%myType=='fdot') then                                                                 ! f_aimDot is prescribed
+     f_aimDot = F_BC%maskFloat * F_BC%values
+   elseif(F_BC%myType=='f') then                                                                    ! aim at end of load case is prescribed
+     f_aimDot = F_BC%maskFloat * (F_BC%values -F_aim)/loadCaseTime
+   endif
+   if (guess) f_aimDot  = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old
+   F_aim_lastInc = F_aim
+
+!--------------------------------------------------------------------------------------------------
+! update coordinates and rate and forward last inc
+   mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(&
+                                            F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)])
+   Fdot =  Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
+                  timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid(3)]))
+   F_lambdaDot =  Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
+                  timeinc_old,guess,F_lambda_lastInc,reshape(F_lambda,[3,3,grid(1),grid(2),grid(3)]))  
+   F_lastInc2 = F_lastInc                
+   F_lastInc     = reshape(F,       [3,3,grid(1),grid(2),grid(3)])
+   F_lambda_lastInc = reshape(F_lambda,[3,3,grid(1),grid(2),grid(3)])
+ endif
+
+ F_aim = F_aim + f_aimDot * timeinc
+ F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, &                                   ! ensure that it matches rotated F_aim
+                                    math_rotate_backward33(F_aim,rotation_BC)), &
+             [9,grid(1),grid(2),grid(3)])
+ F_lambda = reshape(Utilities_forwardField(timeinc,F_lambda_lastInc,F_lambdadot), &
+                    [9,grid(1),grid(2),grid(3)])                                                 ! does not have any average value as boundary condition
+ if (.not. guess) then                                                                              ! large strain forwarding
+   do k = 1_pInt, grid(3); do j = 1_pInt, grid(2); do i = 1_pInt, grid(1)
+      F_lambda33 = reshape(F_lambda(:,i,j,k),[3,3])
+      F_lambda33 = math_mul3333xx33(S_scale,math_mul33x33(F_lambda33, &
+                                  math_mul3333xx33(C_scale,&
+                                                   math_mul33x33(math_transpose33(F_lambda33),&
+                                                                 F_lambda33) -math_I3))*0.5_pReal)&
+                              + math_I3
+      F_lambda(:,i,j,k) = reshape(F_lambda33,[9])
+   enddo; enddo; enddo
+ endif
+ call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr)
+
+end subroutine AL_forward
+
 !--------------------------------------------------------------------------------------------------
 !> @brief destroy routine
 !--------------------------------------------------------------------------------------------------

code/DAMASK_spectral_solverBasicPETSc.f90

@@ -86,7 +86,8 @@ module DAMASK_spectral_SolverBasicPETSc
 
  public :: &
    basicPETSc_init, &
-   basicPETSc_solution ,&
+   basicPETSc_solution, &
+   BasicPETSc_forward, &
    basicPETSc_destroy
  external :: &
    VecDestroy, &
@@ -249,23 +250,13 @@ type(tSolutionState) function basicPETSc_solution( &
  use numerics, only: &
    update_gamma, &
    itmax
- use math, only: &
-   math_mul33x33 ,&
-   math_rotate_backward33
- use mesh, only: &
-   mesh_ipCoordinates,&
-   mesh_deformedCoordsFFT
  use IO, only: &
    IO_write_JobRealFile
  use DAMASK_spectral_Utilities, only: &
    grid, &
-   geomSize, &
    tBoundaryCondition, &
-   Utilities_forwardField, &
-   Utilities_calculateRate, &
    Utilities_maskedCompliance, &
-   Utilities_updateGamma, &
+   Utilities_updateGamma
-   cutBack
  use FEsolving, only: &
    restartWrite, &
    terminallyIll
@@ -282,14 +273,14 @@ type(tSolutionState) function basicPETSc_solution( &
    loadCaseTime, &                                                                                  !< remaining time of current load case
    temperature_bc, &
    density
- logical, intent(in) :: &
-   guess
  type(tBoundaryCondition),      intent(in) :: &
    P_BC, &
    F_BC
  character(len=*), intent(in) :: &
    incInfoIn
  real(pReal), dimension(3,3), intent(in) :: rotation_BC
+ logical, intent(in) :: &
+   guess
  
 !--------------------------------------------------------------------------------------------------
 ! PETSc Data
@@ -298,9 +289,9 @@ type(tSolutionState) function basicPETSc_solution( &
  SNESConvergedReason :: reason
  incInfo = incInfoIn
 
- call DMDAVecGetArrayF90(da,solution_vec,F,ierr)
 !--------------------------------------------------------------------------------------------------
 ! restart information for spectral solver
+ call DMDAVecGetArrayF90(da,solution_vec,F,ierr)
  if (restartWrite) then
    write(6,'(/,a)') ' writing converged results for restart'
    flush(6)
@@ -323,42 +314,9 @@ type(tSolutionState) function basicPETSc_solution( &
    write (777,rec=1) C_volAvgLastInc
    close(777)
  endif 
- mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(&
-                                             F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)])
- if (cutBack) then 
-   F_aim = F_aim_lastInc
-   F = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)]) 
-   C_volAvg = C_volAvgLastInc
- else
-   C_volAvgLastInc = C_volAvg
-
-
-!--------------------------------------------------------------------------------------------------
-! calculate rate for aim
-   if (F_BC%myType=='l') then                                                                       ! calculate f_aimDot from given L and current F
-     f_aimDot = F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim)
-   elseif(F_BC%myType=='fdot') then                                                                 ! f_aimDot is prescribed
-     f_aimDot = F_BC%maskFloat * F_BC%values
-   elseif(F_BC%myType=='f') then                                                                    ! aim at end of load case is prescribed
-     f_aimDot = F_BC%maskFloat * (F_BC%values -F_aim)/loadCaseTime
-  endif
-   if (guess) f_aimDot  = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old
-   F_aim_lastInc = F_aim
-
-!--------------------------------------------------------------------------------------------------
-! update coordinates and rate and forward last inc
-   mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(&
-                                            F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)])
-   Fdot =  Utilities_calculateRate(math_rotate_backward33(f_aimDot,params%rotation_BC), &
-                                 timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid(3)]))
-   F_lastInc2 = F_lastInc
-   F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid(3)])
- endif
- F_aim = F_aim + f_aimDot * timeinc
-
- F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot,math_rotate_backward33(F_aim, &
-                                                            rotation_BC)),[9,grid(1),grid(2),grid(3)])
  call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
+ 
+ BasicPETSc_solution%converged =.false.
   
 !--------------------------------------------------------------------------------------------------
 ! update stiffness (and gamma operator)
@@ -607,6 +565,78 @@ subroutine BasicPETSc_convergedKSP(ksp_local,PETScIter,fnorm,reason,dummy,ierr)
  
 end subroutine BasicPETSc_convergedKSP
 
+!--------------------------------------------------------------------------------------------------
+!> @brief forwarding routine
+!--------------------------------------------------------------------------------------------------
+subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_BC)
+ use math, only: &
+   math_mul33x33 ,&
+   math_rotate_backward33
+ use DAMASK_spectral_Utilities, only: &
+   grid, &
+   geomSize, &
+   Utilities_calculateRate, &
+   Utilities_forwardField, &
+   tBoundaryCondition, &
+   cutBack
+ use mesh, only: &
+   mesh_ipCoordinates,&
+   mesh_deformedCoordsFFT
+
+ implicit none
+#include <finclude/petscdmda.h90>
+ real(pReal), intent(in) :: &
+   timeinc_old, &
+   timeinc, &
+   loadCaseTime                                                                                     !< remaining time of current load case
+ type(tBoundaryCondition),      intent(in) :: &
+   P_BC, &
+   F_BC
+ real(pReal), dimension(3,3), intent(in) :: rotation_BC
+ logical, intent(in) :: &
+   guess
+ PetscScalar, pointer :: F(:,:,:,:)
+ PetscErrorCode :: ierr
+
+!--------------------------------------------------------------------------------------------------
+! update coordinates and rate and forward last inc
+ call DMDAVecGetArrayF90(da,solution_vec,F,ierr)
+ if (cutBack) then 
+   F_aim = F_aim_lastInc
+   F    = reshape(F_lastInc,    [9,grid(1),grid(2),grid(3)]) 
+   C_volAvg = C_volAvgLastInc
+ else
+   C_volAvgLastInc = C_volAvg
+!--------------------------------------------------------------------------------------------------
+! calculate rate for aim
+   if (F_BC%myType=='l') then                                                                       ! calculate f_aimDot from given L and current F
+     f_aimDot = F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim)
+   elseif(F_BC%myType=='fdot') then                                                                 ! f_aimDot is prescribed
+     f_aimDot = F_BC%maskFloat * F_BC%values
+   elseif(F_BC%myType=='f') then                                                                    ! aim at end of load case is prescribed
+     f_aimDot = F_BC%maskFloat * (F_BC%values -F_aim)/loadCaseTime
+   endif
+   if (guess) f_aimDot  = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old
+   F_aim_lastInc = F_aim
+
+!--------------------------------------------------------------------------------------------------
+! update coordinates and rate and forward last inc
+   mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(&
+                                            F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)])
+   Fdot =  Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
+                  timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid(3)]))
+   F_lastInc2 = F_lastInc
+   F_lastInc     = reshape(F,       [3,3,grid(1),grid(2),grid(3)])
+ endif
+ F_aim = F_aim + f_aimDot * timeinc
+
+!--------------------------------------------------------------------------------------------------
+! update local deformation gradient
+ F     = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, &                                      ! ensure that it matches rotated F_aim
+                               math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid(3)])
+ call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
+
+end subroutine BasicPETSc_forward
 
 !--------------------------------------------------------------------------------------------------
 !> @brief destroy routine

code/DAMASK_spectral_solverPolarisation.f90

@@ -100,6 +100,7 @@ module DAMASK_spectral_solverPolarisation
  public :: &
    Polarisation_init, &
    Polarisation_solution, &
+   Polarisation_forward, &
    Polarisation_destroy
  external :: &
    VecDestroy, &
@@ -277,25 +278,15 @@ type(tSolutionState) function &
  use numerics, only: &
    update_gamma
  use math, only: &
-   math_mul33x33 ,&
+   math_invSym3333
-   math_mul3333xx33, &
-   math_rotate_backward33, &
-   math_invSym3333, &
-   math_transpose33
-use mesh, only: &
-   mesh_ipCoordinates, &
-   mesh_deformedCoordsFFT
  use IO, only: &
    IO_write_jobRealFile
  use DAMASK_spectral_Utilities, only: &
    grid, &
    geomSize, &
    tBoundaryCondition, &
-   Utilities_forwardField, &
-   Utilities_calculateRate, &
    Utilities_maskedCompliance, &
-   Utilities_updateGamma, &
+   Utilities_updateGamma
-   cutBack
  use FEsolving, only: &
    restartWrite, &
    terminallyIll
@@ -327,16 +318,13 @@ use mesh, only: &
  PetscErrorCode :: ierr   
  SNESConvergedReason :: reason
 
- integer(pInt) :: i, j, k
- real(pReal), dimension(3,3) :: F_lambda33
-
  incInfo = incInfoIn
- call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr)
- F => xx_psc(0:8,:,:,:)
- F_tau => xx_psc(9:17,:,:,:)
  
 !--------------------------------------------------------------------------------------------------
 ! restart information for spectral solver
+ call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr)
+ F => xx_psc(0:8,:,:,:)
+ F_tau => xx_psc(9:17,:,:,:)
  if (restartWrite) then
    write(6,'(/,a)') ' writing converged results for restart'
    flush(6)
@@ -365,60 +353,9 @@ use mesh, only: &
    write (777,rec=1) C_volAvgLastInc
    close(777)
  endif 
+ call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr)
  Polarisation_solution%converged =.false.
 
- if (cutBack) then 
-   F_aim = F_aim_lastInc
-   F_tau= reshape(F_tau_lastInc,[9,grid(1),grid(2),grid(3)]) 
-   F    = reshape(F_lastInc,    [9,grid(1),grid(2),grid(3)]) 
-   C_volAvg = C_volAvgLastInc
- else
-   C_volAvgLastInc = C_volAvg
-!--------------------------------------------------------------------------------------------------
-! calculate rate for aim
-   if (F_BC%myType=='l') then                                                                       ! calculate f_aimDot from given L and current F
-     f_aimDot = F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim)
-   elseif(F_BC%myType=='fdot') then                                                                 ! f_aimDot is prescribed
-     f_aimDot = F_BC%maskFloat * F_BC%values
-   elseif(F_BC%myType=='f') then                                                                    ! aim at end of load case is prescribed
-     f_aimDot = F_BC%maskFloat * (F_BC%values -F_aim)/loadCaseTime
-   endif
-   if (guess) f_aimDot  = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old
-   F_aim_lastInc = F_aim
-
-!--------------------------------------------------------------------------------------------------
-! update coordinates and rate and forward last inc
-   mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(&
-                                            F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)])
-   Fdot =  Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
-                  timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid(3)]))
-   F_tauDot =  Utilities_calculateRate(math_rotate_backward33(2.0_pReal*f_aimDot,rotation_BC), &
-                  timeinc_old,guess,F_tau_lastInc,reshape(F_tau,[3,3,grid(1),grid(2),grid(3)]))  
-   F_lastInc2 = F_lastInc
-   F_lastInc     = reshape(F,       [3,3,grid(1),grid(2),grid(3)])
-   F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid(3)])
- endif
- F_aim = F_aim + f_aimDot * timeinc
-
-!--------------------------------------------------------------------------------------------------
-! update local deformation gradient
- F     = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, &                                      ! ensure that it matches rotated F_aim
-                               math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid(3)])
- F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot),  [9,grid(1),grid(2),grid(3)]) ! does not have any average value as boundary condition
- if (.not. guess) then                                                                                 ! large strain forwarding
-   do k = 1_pInt, grid(3); do j = 1_pInt, grid(2); do i = 1_pInt, grid(1)
-      F_lambda33 = reshape(F_tau(:,i,j,k)-F(:,i,j,k),[3,3])
-      F_lambda33 = math_mul3333xx33(S_scale,math_mul33x33(F_lambda33, &
-                                  math_mul3333xx33(C_scale,&
-                                                   math_mul33x33(math_transpose33(F_lambda33),&
-                                                                 F_lambda33) -math_I3))*0.5_pReal)&
-                              + math_I3
-      F_tau(:,i,j,k) = reshape(F_lambda33,[9])+F(:,i,j,k)
-   enddo; enddo; enddo
- endif
- call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr)
- CHKERRQ(ierr)
-
 !--------------------------------------------------------------------------------------------------
 ! update stiffness (and gamma operator)
  S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C_volAvg)
@@ -694,6 +631,101 @@ subroutine Polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,
 
 end subroutine Polarisation_converged
 
+!--------------------------------------------------------------------------------------------------
+!> @brief forwarding routine
+!--------------------------------------------------------------------------------------------------
+subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_BC)
+ use math, only: &
+   math_mul33x33, &
+   math_mul3333xx33, &
+   math_transpose33, &
+   math_rotate_backward33
+ use DAMASK_spectral_Utilities, only: &
+   grid, &
+   geomSize, &
+   Utilities_calculateRate, &
+   Utilities_forwardField, &
+   tBoundaryCondition, &
+   cutBack
+ use mesh, only: &
+   mesh_ipCoordinates,&
+   mesh_deformedCoordsFFT
+
+ implicit none
+#include <finclude/petscdmda.h90>
+ real(pReal), intent(in) :: &
+   timeinc_old, &
+   timeinc, &
+   loadCaseTime                                                                                     !< remaining time of current load case
+ type(tBoundaryCondition),      intent(in) :: &
+   P_BC, &
+   F_BC
+ real(pReal), dimension(3,3), intent(in) :: rotation_BC
+ logical, intent(in) :: &
+   guess
+ PetscErrorCode :: ierr
+ PetscScalar, dimension(:,:,:,:), pointer :: xx_psc, F, F_tau
+ integer(pInt) :: i, j, k
+ real(pReal), dimension(3,3) :: F_lambda33
+
+!--------------------------------------------------------------------------------------------------
+! update coordinates and rate and forward last inc
+ call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr)
+ F => xx_psc(0:8,:,:,:)
+ F_tau => xx_psc(9:17,:,:,:)
+ if (cutBack) then 
+   F_aim = F_aim_lastInc
+   F_tau= reshape(F_tau_lastInc,[9,grid(1),grid(2),grid(3)]) 
+   F    = reshape(F_lastInc,    [9,grid(1),grid(2),grid(3)]) 
+   C_volAvg = C_volAvgLastInc
+ else
+   C_volAvgLastInc = C_volAvg
+!--------------------------------------------------------------------------------------------------
+! calculate rate for aim
+   if (F_BC%myType=='l') then                                                                       ! calculate f_aimDot from given L and current F
+     f_aimDot = F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim)
+   elseif(F_BC%myType=='fdot') then                                                                 ! f_aimDot is prescribed
+     f_aimDot = F_BC%maskFloat * F_BC%values
+   elseif(F_BC%myType=='f') then                                                                    ! aim at end of load case is prescribed
+     f_aimDot = F_BC%maskFloat * (F_BC%values -F_aim)/loadCaseTime
+   endif
+   if (guess) f_aimDot  = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old
+   F_aim_lastInc = F_aim
+
+!--------------------------------------------------------------------------------------------------
+! update coordinates and rate and forward last inc
+   mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(&
+                                            F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)])
+   Fdot =  Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
+                  timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid(3)]))
+   F_tauDot =  Utilities_calculateRate(math_rotate_backward33(2.0_pReal*f_aimDot,rotation_BC), &
+                  timeinc_old,guess,F_tau_lastInc,reshape(F_tau,[3,3,grid(1),grid(2),grid(3)]))  
+   F_lastInc2 = F_lastInc
+   F_lastInc     = reshape(F,       [3,3,grid(1),grid(2),grid(3)])
+   F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid(3)])
+ endif
+ F_aim = F_aim + f_aimDot * timeinc
+
+!--------------------------------------------------------------------------------------------------
+! update local deformation gradient
+ F     = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, &                                      ! ensure that it matches rotated F_aim
+                               math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid(3)])
+ F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot),  [9,grid(1),grid(2),grid(3)]) ! does not have any average value as boundary condition
+ if (.not. guess) then                                                                                 ! large strain forwarding
+   do k = 1_pInt, grid(3); do j = 1_pInt, grid(2); do i = 1_pInt, grid(1)
+      F_lambda33 = reshape(F_tau(:,i,j,k)-F(:,i,j,k),[3,3])
+      F_lambda33 = math_mul3333xx33(S_scale,math_mul33x33(F_lambda33, &
+                                  math_mul3333xx33(C_scale,&
+                                                   math_mul33x33(math_transpose33(F_lambda33),&
+                                                                 F_lambda33) -math_I3))*0.5_pReal)&
+                              + math_I3
+      F_tau(:,i,j,k) = reshape(F_lambda33,[9])+F(:,i,j,k)
+   enddo; enddo; enddo
+ endif
+ call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr)
+
+end subroutine Polarisation_forward
+
 !--------------------------------------------------------------------------------------------------
 !> @brief destroy routine
 !--------------------------------------------------------------------------------------------------