restart is working for all solvers, but slight deviations still occur. Reason is most probably, that for CPFEM init, a 0s timestep is required which seems to alter the state a tiny bit, leading to small deviations. For AL and Polarization, this is even more severe since the regular call to CPFEM has F-residual as the current deformation gradient, however the init step uses only F.

code/DAMASK_spectral_driver.f90

@@ -355,18 +355,14 @@ program DAMASK_spectral_Driver
  if (appendToOutFile) then                                                                          ! after restart, append to existing results file
    open(newunit=resUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//&
                                '.spectralOut',form='UNFORMATTED', position='APPEND', status='OLD')
-                             !   '.spectralOut',access='STREAM', position='APPEND', status='OLD')
    open(newunit=statUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//&
                                '.sta',form='FORMATTED', position='APPEND', status='OLD')
  else                                                                                               ! open new files ...
    open(newunit=resUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//&
                                '.spectralOut',form='UNFORMATTED',status='REPLACE')
-                             !   '.spectralOut',access='STREAM',status='REPLACE')
    write(resUnit) 'load',       trim(loadCaseFile)                                                  ! ... and write header
    write(resUnit) 'workingdir', trim(getSolverWorkingDirectoryName())
    write(resUnit) 'geometry',   trim(geometryFile)
-   !write(resUnit) 'grid', grid
-   !write(resUnit) 'size',  geomSize
    write(resUnit) 'resolution', grid
    write(resUnit) 'dimension',  geomSize
    write(resUnit) 'materialpoint_sizeResults', materialpoint_sizeResults

code/DAMASK_spectral_solverAL.f90

@@ -163,9 +163,6 @@ subroutine AL_init(temperature)
  real(pReal), dimension(:,:,:,:,:), allocatable :: P
  real(pReal), dimension(3,3) :: &
    temp33_Real = 0.0_pReal
- real(pReal), dimension(3,3,3,3) :: &
-   temp3333_Real = 0.0_pReal, &
-   temp3333_Real2 = 0.0_pReal
 
  PetscErrorCode :: ierr
  PetscObject :: dummy
@@ -214,7 +211,8 @@ subroutine AL_init(temperature)
    F = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)])
    F_lambda = F
    F_lambda_lastInc = F_lastInc
- elseif (restartInc > 1_pInt) then 
+
+ elseif (restartInc > 1_pInt) then                                                                  ! read in F to calculate coordinates and initialize CPFEM general
    if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) &
      write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
      'reading values of increment ', restartInc - 1_pInt, ' from file'
@@ -244,28 +242,35 @@ subroutine AL_init(temperature)
    call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(f_aimDot))
    read (777,rec=1) f_aimDot
    close (777)
+ endif
+ 
+ mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(&
+                              F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)])
+ call Utilities_constitutiveResponse(F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid(3)]),&
+   temperature,0.0_pReal,P,C_volAvg,C_minMaxAvg,temp33_Real,.false.,math_I3)
+ nullify(F)
+ nullify(F_lambda)
+ call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr)                            ! write data back to PETSc
+                             
+ if (restartInc > 1_pInt) then                                                                      ! using old values from files
+   if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) &
+     write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
+     'reading more values of increment', restartInc - 1_pInt, 'from file'
+   flush(6)
+   call IO_read_realFile(777,'F_lambda',trim(getSolverJobName()),size(F_lambda))
+   read (777,rec=1) F_lambda
    call IO_read_realFile(777,'C_volAvg',trim(getSolverJobName()),size(C_volAvg))
    read (777,rec=1) C_volAvg
    close (777)
    call IO_read_realFile(777,'C_volAvgLastInc',trim(getSolverJobName()),size(C_volAvgLastInc))
    read (777,rec=1) C_volAvgLastInc
    close (777)
-   call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(temp3333_Real))
+   call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg))
    read (777,rec=1) C_minMaxAvg
    close (777)
  endif
- mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(&
-                                              F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)])
- call Utilities_constitutiveResponse(F,F,temperature,0.0_pReal,P,temp3333_Real,temp3333_Real2,&
-                                temp33_Real,.false.,math_I3)
- call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr)
 
-!--------------------------------------------------------------------------------------------------
+
-! reference stiffness
- if (restartInc == 1_pInt) then                                                                     ! use initial stiffness as reference stiffness
-   C_minMaxAvg = temp3333_Real2
-   C_volAvg = temp3333_Real
- endif 
 
  call Utilities_updateGamma(C_minMaxAvg,.True.)
  C_scale = C_minMaxAvg
@@ -419,11 +424,11 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
  integer(pInt) :: &
    i, j, k, e
 
- F              => x_scal(1:3,1:3,1,&
+ F                => x_scal(1:3,1:3,1,&
   XG_RANGE,YG_RANGE,ZG_RANGE)
- F_lambda          => x_scal(1:3,1:3,2,&
+ F_lambda         => x_scal(1:3,1:3,2,&
   XG_RANGE,YG_RANGE,ZG_RANGE)
- residual_F     => f_scal(1:3,1:3,1,&
+ residual_F       => f_scal(1:3,1:3,1,&
   X_RANGE,Y_RANGE,Z_RANGE)
  residual_F_lambda => f_scal(1:3,1:3,2,&
   X_RANGE,Y_RANGE,Z_RANGE)
@@ -434,7 +439,7 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
  F_av = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt
  
  if(nfuncs== 0 .and. PETScIter == 0) totalIter = -1_pInt                                            ! new increment
- if (totalIter <= PETScIter) then                                                                   ! new iteration
+ if(totalIter <= PETScIter) then                                                                    ! new iteration
 !--------------------------------------------------------------------------------------------------
 ! report begin of new iteration
    totalIter = totalIter + 1_pInt
@@ -507,7 +512,7 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
  e = 0_pInt
  do k = 1_pInt, grid(3); 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,e) + C_scale), &
+   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_lambda(1:3,1:3,i,j,k) - math_I3))) &
@@ -675,6 +680,7 @@ subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_
    call IO_write_jobRealFile(777,'C_volAvgLastInc',size(C_volAvgLastInc))
    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)])
@@ -682,13 +688,13 @@ subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_
  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)]) 
+   F        = reshape(F_lastInc,       [9,grid(1),grid(2),grid(3)]) 
    C_volAvg = C_volAvgLastInc
  else
    ForwardData = .True.
+   C_volAvgLastInc = C_volAvg
 !--------------------------------------------------------------------------------------------------
 ! 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
@@ -707,8 +713,8 @@ subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_
                   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_lastInc2       = F_lastInc                
-   F_lastInc     = reshape(F,       [3,3,grid(1),grid(2),grid(3)])
+   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
 
@@ -720,13 +726,13 @@ subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_
                     [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 = reshape(F_lambda(1:9,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])
+      F_lambda(1:9,i,j,k) = reshape(F_lambda33,[9])
    enddo; enddo; enddo
  endif
  call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr)

code/DAMASK_spectral_solverBasic.f90

@@ -119,7 +119,7 @@ subroutine basic_init(temperature)
  elseif (restartInc > 1_pInt) then                                                                  ! using old values from file                                                      
    if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) &
      write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
-     'reading deformation gradients of increment', restartInc - 1_pInt, 'from file'
+     'reading values of increment ', restartInc - 1_pInt, ' from file'
    flush(6)
    call IO_read_realFile(777,'F',trim(getSolverJobName()),size(F))
    read (777,rec=1) F
@@ -127,31 +127,30 @@ subroutine basic_init(temperature)
    call IO_read_realFile(777,'F_lastInc',trim(getSolverJobName()),size(F_lastInc))
    read (777,rec=1) F_lastInc
    close (777)
+   call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(f_aimDot))
+   read (777,rec=1) f_aimDot
+   close (777)
+   F_aim         = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt                                         ! average of F
+   F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt                                 ! average of F_lastInc 
  endif
 
  mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,F),[3,1,product(grid)])
- call Utilities_constitutiveResponse(F,F,temperature,0.0_pReal,P,C,C_minmaxAvg,&
+ call Utilities_constitutiveResponse(F_lastInc,F,temperature,0.0_pReal,P,C,C_minMaxAvg,&
                                      temp33_Real,.false.,math_I3)                                   ! constitutive response with no deformation in no time to get reference stiffness
    
- if (restartInc > 1_pInt) then                                                                      ! using old values from file                                                      
+ if (restartInc > 1_pInt) then                                                                      ! using old values from files
-   F_aim         = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt                                         ! average of F
-   F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt                                 ! average of F_lastInc 
-
    if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) &
      write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
      'reading more values of increment', restartInc - 1_pInt, 'from file'
    flush(6)
-   call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(f_aimDot))
-   read (777,rec=1) f_aimDot
-   close (777)
    call IO_read_realFile(777,'C',trim(getSolverJobName()),size(C))
    read (777,rec=1) C
    close (777)
    call IO_read_realFile(777,'C_lastInc',trim(getSolverJobName()),size(C_lastInc))
    read (777,rec=1) C_lastInc
    close (777)
-   call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minmaxAvg))
+   call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg))
-   read (777,rec=1) C_minmaxAvg
+   read (777,rec=1) C_minMaxAvg
    close (777)
  endif
    

code/DAMASK_spectral_solverBasicPETSc.f90

@@ -198,7 +198,7 @@ subroutine basicPETSc_init(temperature)
  elseif (restartInc > 1_pInt) then                                                                  ! using old values from file                                                      
    if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) &
      write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
-     'reading deformation gradients of increment', restartInc - 1_pInt, 'from file'
+     'reading values of increment ', restartInc - 1_pInt, ' from file'
    flush(6)
    call IO_read_realFile(777,'F',trim(getSolverJobName()),size(F))
    read (777,rec=1) F
@@ -209,35 +209,33 @@ subroutine basicPETSc_init(temperature)
    call IO_read_realFile(777,'F_lastInc2',trim(getSolverJobName()),size(F_lastInc2))
    read (777,rec=1) F_lastInc2
    close (777)
+   call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(f_aimDot))
+   read (777,rec=1) f_aimDot
+   close (777)
+   F_aim         = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3])                         ! average of F
+   F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt                                 ! average of F_lastInc 
  endif
 
  mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(&
                                               F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)])
- call Utilities_constitutiveResponse(&
+ call Utilities_constitutiveResponse(F_lastInc, &
     reshape(F(0:8,0:grid(1)-1_pInt,0:grid(2)-1_pInt,0:grid(3)-1_pInt),[3,3,grid(1),grid(2),grid(3)]),&
-    reshape(F(0:8,0:grid(1)-1_pInt,0:grid(2)-1_pInt,0:grid(3)-1_pInt),[3,3,grid(1),grid(2),grid(3)]),&
+    temperature,0.0_pReal,P,C_volAvg,C_minMaxAvg,temp33_Real,.false.,math_I3)
-    temperature,0.0_pReal,P,C_volAvg,C_minmaxAvg,temp33_Real,.false.,math_I3)
+ call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)                                 ! write data back to PETSc
- call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)                                 ! write data back into PETSc
-
- if (restartInc > 1_pInt) then                                                                      ! using old values from file                                                      
-   F_aim         = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3])                         ! average of F
-   F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt                                 ! average of F_lastInc 
 
+ if (restartInc > 1_pInt) then                                                                      ! using old values from files                                                    
    if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) &
      write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
      'reading more values of increment', restartInc - 1_pInt, 'from file'
    flush(6)
-   call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(f_aimDot))
-   read (777,rec=1) f_aimDot
-   close (777)
    call IO_read_realFile(777,'C_volAvg',trim(getSolverJobName()),size(C_volAvg))
    read (777,rec=1) C_volAvg
    close (777)
    call IO_read_realFile(777,'C_volAvgLastInc',trim(getSolverJobName()),size(C_volAvgLastInc))
    read (777,rec=1) C_volAvgLastInc
    close (777)
-   call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minmaxAvg))
+   call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg))
-   read (777,rec=1) C_minmaxAvg
+   read (777,rec=1) C_minMaxAvg
    close (777)
  endif 
    

code/DAMASK_spectral_solverPolarisation.f90

@@ -163,9 +163,6 @@ subroutine Polarisation_init(temperature)
  real(pReal), dimension(:,:,:,:,:), allocatable :: P
  real(pReal), dimension(3,3) :: &
    temp33_Real = 0.0_pReal
- real(pReal), dimension(3,3,3,3) :: &
-   temp3333_Real = 0.0_pReal, &
-   temp3333_Real2 = 0.0_pReal
 
  PetscErrorCode :: ierr
  PetscObject :: dummy
@@ -211,9 +208,9 @@ subroutine Polarisation_init(temperature)
  if (restartInc == 1_pInt) then                                                                     ! no deformation (no restart)
    F_lastInc     = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid(3))                    ! initialize to identity
    F_lastInc2 = F_lastInc
-   F_tau_lastInc = 2.0_pReal*F_lastInc
    F = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)])
    F_tau = 2.0_pReal* F
+   F_tau_lastInc = 2.0_pReal*F_lastInc
  elseif (restartInc > 1_pInt) then 
    if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) &
      write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
@@ -228,8 +225,6 @@ subroutine Polarisation_init(temperature)
    call IO_read_realFile(777,'F_lastInc2',trim(getSolverJobName()),size(F_lastInc2))
    read (777,rec=1) F_lastInc2
    close (777)
-   F_aim         = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3])                         ! average of F
-   F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt                                 ! average of F_lastInc 
    call IO_read_realFile(777,'F_tau',trim(getSolverJobName()),size(F_tau))
    read (777,rec=1) F_tau
    close (777)
@@ -246,32 +241,35 @@ subroutine Polarisation_init(temperature)
    call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(f_aimDot))
    read (777,rec=1) f_aimDot
    close (777)
+ endif
+
+ mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(&
+                              F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)])
+ call Utilities_constitutiveResponse(F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid(3)]),&
+   temperature,0.0_pReal,P,C_volAvg,C_minMaxAvg,temp33_Real,.false.,math_I3)
+ nullify(F)
+ nullify(F_tau)
+ call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr)                            ! write data back to PETSc
+
+ if (restartInc > 1_pInt) then                                                                      ! using old values from files
+   if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) &
+     write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
+     'reading more values of increment', restartInc - 1_pInt, 'from file'
+   flush(6)
    call IO_read_realFile(777,'C_volAvg',trim(getSolverJobName()),size(C_volAvg))
    read (777,rec=1) C_volAvg
    close (777)
    call IO_read_realFile(777,'C_volAvgLastInc',trim(getSolverJobName()),size(C_volAvgLastInc))
    read (777,rec=1) C_volAvgLastInc
    close (777)
-   call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(temp3333_Real))
+   call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg))
    read (777,rec=1) C_minMaxAvg
    close (777)
  endif
- mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(&
-                                              F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)])
- call Utilities_constitutiveResponse(F,F,temperature,0.0_pReal,P,temp3333_Real,temp3333_Real2,&
-                                temp33_Real,.false.,math_I3)
- call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr)
 
-!--------------------------------------------------------------------------------------------------
+ call Utilities_updateGamma(C_minMaxAvg,.True.)
-! reference stiffness
+ C_scale = C_minMaxAvg
- if (restartInc == 1_pInt) then                                                                     ! use initial stiffness as reference stiffness
+ S_scale = math_invSym3333(C_minMaxAvg)
-   C_minMaxAvg = temp3333_Real2
-   C_volAvg = temp3333_Real
- endif 
-
- call Utilities_updateGamma(temp3333_Real2,.True.)
- C_scale = temp3333_Real2
- S_scale = math_invSym3333(temp3333_Real2)
  
 end subroutine Polarisation_init
 
@@ -508,7 +506,7 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
  do k = 1_pInt, grid(3); 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,e) + C_scale), &
+     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)
@@ -720,13 +718,13 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC
  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 = reshape(F_tau(1:9,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)
+      F_tau(1:9,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)

code/DAMASK_spectral_utilities.f90

@@ -862,6 +862,7 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,&
  if (forwardData) then                                                                              ! aging results
    calcMode    = ior(calcMode,    CPFEM_AGERESULTS)             
    collectMode = ior(collectMode, CPFEM_BACKUPJACOBIAN)
+   materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid)])
  endif
  if (cutBack) then                                                                                  ! restore saved variables
   collectMode = ior(collectMode , CPFEM_RESTOREJACOBIAN)
@@ -872,8 +873,7 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,&
  call CPFEM_general(collectMode,usePingPong,F_lastInc(1:3,1:3,1,1,1),F(1:3,1:3,1,1,1), &            ! collect mode handles Jacobian backup / restoration
                    crystallite_temperature(1,1),timeinc,1_pInt,1_pInt)
  
- materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid)])
+ materialpoint_F  = reshape(F,[3,3,1,product(grid)])
- materialpoint_F  = reshape(F,         [3,3,1,product(grid)])
 
  call debug_reset()