changed up handling of delta states, and some bug fixes

code/DAMASK_spectral_solverAL.f90

@@ -13,7 +13,8 @@ module DAMASK_spectral_solverAL
  use math, only: &
    math_I3
  use DAMASK_spectral_utilities, only: &
-   tSolutionState
+   tSolutionState, &
+   tSolutionParams
  
  implicit none
  private
@@ -24,14 +25,6 @@ module DAMASK_spectral_solverAL
    
 !--------------------------------------------------------------------------------------------------
 ! derived types 
- type tSolutionParams                                                                               !< @todo use here the type definition for a full loadcase including mask
-   real(pReal), dimension(3,3) :: P_BC, rotation_BC
-   real(pReal) :: timeinc
-   real(pReal) :: timeincOld
-   real(pReal) :: temperature
-   real(pReal) :: density
- end type tSolutionParams
- 
  type(tSolutionParams), private :: params
  real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal
 
@@ -96,9 +89,7 @@ module DAMASK_spectral_solverAL
    SNESSetConvergenceTest, &
    SNESSetFromOptions, &
    SNESCreate, &
-   MPI_Abort, &
+   MPI_Abort
-   MPI_Bcast, &
-   MPI_Allreduce
 
 contains
 
@@ -124,12 +115,9 @@ subroutine AL_init(temperature)
  use DAMASK_interface, only: &
    getSolverJobName
  use DAMASK_spectral_Utilities, only: &
-   Utilities_init, &
    Utilities_constitutiveResponse, &
    Utilities_updateGamma, &
-   Utilities_updateIPcoords, &
+   Utilities_updateIPcoords
-   grid1Red, &
-   wgt
  use mesh, only: &
    gridLocal, &
    gridGlobal
@@ -150,7 +138,6 @@ subroutine AL_init(temperature)
  integer(pInt) :: proc
  character(len=1024) :: rankStr
  
- call Utilities_init()
  if (worldrank == 0_pInt) then
    write(6,'(/,a)') ' <<<+-  DAMASK_spectral_solverAL init  -+>>>'
    write(6,'(a)') ' $Id$'
@@ -169,6 +156,7 @@ subroutine AL_init(temperature)
 !--------------------------------------------------------------------------------------------------
 ! PETSc Init
  call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr)
+ call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) 
  allocate(localK(worldsize), source = 0); localK(worldrank+1) = gridLocal(3)
  do proc = 1, worldsize
    call MPI_Bcast(localK(proc),1,MPI_INTEGER,proc-1,PETSC_COMM_WORLD,ierr)
@@ -182,6 +170,7 @@ subroutine AL_init(temperature)
         gridLocal (1),gridLocal (2),localK, &                                                       ! local grid
         da,ierr)                                                                                    ! handle, error
  CHKERRQ(ierr)
+ call SNESSetDM(snes,da,ierr); CHKERRQ(ierr)
  call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr)
  call DMDASNESSetFunctionLocal(da,INSERT_VALUES,AL_formResidual,dummy,ierr)
  CHKERRQ(ierr)
@@ -266,7 +255,7 @@ end subroutine AL_init
 !--------------------------------------------------------------------------------------------------
 type(tSolutionState) function &
   AL_solution(incInfoIn,guess,timeinc,timeinc_old,loadCaseTime,P_BC,F_BC,temperature_bc, &
-                                                                     rotation_BC,density)
+                                                                     rotation_BC)
  use numerics, only: &
    update_gamma
  use math, only: &
@@ -287,8 +276,7 @@ type(tSolutionState) function &
    timeinc, &                                                                                       !< increment in time for current solution
    timeinc_old, &                                                                                   !< increment in time of last increment
    loadCaseTime, &                                                                                  !< remaining time of current load case
-   temperature_bc, &
+   temperature_bc
-   density
  logical, intent(in) :: &
    guess
  type(tBoundaryCondition),      intent(in) :: &
@@ -324,7 +312,6 @@ type(tSolutionState) function &
  params%timeinc = timeinc
  params%timeincOld = timeinc_old
  params%temperature = temperature_bc
- params%density = density
 
 !--------------------------------------------------------------------------------------------------
 ! solve BVP 
@@ -363,16 +350,14 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
    math_transpose33, &
    math_mul3333xx33, &
    math_invSym3333, &
-   math_mul33x33, &
+   math_mul33x33
-   PI
  use DAMASK_spectral_Utilities, only: &
    wgt, &
-   field_realMPI, &
+   tensorField_realMPI, &
-   field_fourierMPI, &
+   utilities_FFTtensorForward, &
-   Utilities_FFTforward, &
+   utilities_fourierGammaConvolution, &
-   Utilities_fourierConvolution, &
    Utilities_inverseLaplace, &
-   Utilities_FFTbackward, &
+   utilities_FFTtensorBackward, &
    Utilities_constitutiveResponse, &
    Utilities_divergenceRMS, &
    Utilities_curlRMS
@@ -444,9 +429,9 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
 
 !--------------------------------------------------------------------------------------------------
 ! 
- field_realMPI = 0.0_pReal
+ tensorField_realMPI = 0.0_pReal
  do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt, gridLocal(1)
-   field_realMPI(1:3,1:3,i,j,k) = &
+   tensorField_realMPI(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_lambda(1:3,1:3,i,j,k) - math_I3))
@@ -455,13 +440,13 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
  
 !--------------------------------------------------------------------------------------------------
 ! doing convolution in Fourier space 
- call Utilities_FFTforward()
+ call utilities_FFTtensorForward()
- call Utilities_fourierConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC)) 
+ call utilities_fourierGammaConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC)) 
- call Utilities_FFTbackward()
+ call utilities_FFTtensorBackward()
 
 !--------------------------------------------------------------------------------------------------
 ! constructing residual                         
- residual_F_lambda = polarBeta*F - field_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3))
+ residual_F_lambda = polarBeta*F - tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3))
 
 !--------------------------------------------------------------------------------------------------
 ! evaluate constitutive response
@@ -473,11 +458,11 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
 
 !--------------------------------------------------------------------------------------------------
 ! calculate divergence
- field_realMPI = 0.0_pReal
+ tensorField_realMPI = 0.0_pReal
- field_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = residual_F
+ tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = residual_F
- call Utilities_FFTforward()
+ call utilities_FFTtensorForward()
  err_div = Utilities_divergenceRMS()
- call Utilities_FFTbackward()
+ call utilities_FFTtensorBackward()
  
 !--------------------------------------------------------------------------------------------------
 ! constructing residual
@@ -493,11 +478,11 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
  
 !--------------------------------------------------------------------------------------------------
 ! calculating curl 
- field_realMPI = 0.0_pReal
+ tensorField_realMPI = 0.0_pReal
- field_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = F
+ tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = F
- call Utilities_FFTforward()
+ call utilities_FFTtensorForward()
  err_curl = Utilities_curlRMS()
- call Utilities_FFTbackward()
+ call utilities_FFTtensorBackward()
  
 end subroutine AL_formResidual
 
@@ -727,7 +712,6 @@ subroutine AL_destroy()
  call VecDestroy(solution_vec,ierr); CHKERRQ(ierr)
  call SNESDestroy(snes,ierr); CHKERRQ(ierr)
  call DMDestroy(da,ierr); CHKERRQ(ierr)
- call Utilities_destroy()
 
 end subroutine AL_destroy
 

code/DAMASK_spectral_solverPolarisation.f90

@@ -13,7 +13,8 @@ module DAMASK_spectral_solverPolarisation
  use math, only: &
    math_I3
  use DAMASK_spectral_utilities, only: &
-   tSolutionState
+   tSolutionState, &
+   tSolutionParams
  
  implicit none
  private
@@ -24,14 +25,6 @@ module DAMASK_spectral_solverPolarisation
    
 !--------------------------------------------------------------------------------------------------
 ! derived types 
- type tSolutionParams                                                                               !< @todo use here the type definition for a full loadcase including mask
-   real(pReal), dimension(3,3) :: P_BC, rotation_BC
-   real(pReal) :: timeinc
-   real(pReal) :: timeincOld
-   real(pReal) :: temperature
-   real(pReal) :: density
- end type tSolutionParams
- 
  type(tSolutionParams), private :: params
  real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal
 
@@ -124,12 +117,9 @@ subroutine Polarisation_init(temperature)
  use DAMASK_interface, only: &
    getSolverJobName
  use DAMASK_spectral_Utilities, only: &
-   Utilities_init, &
    Utilities_constitutiveResponse, &
    Utilities_updateGamma, &
-   Utilities_updateIPcoords, &
+   Utilities_updateIPcoords
-   grid1Red, &
-   wgt
  use mesh, only: &
    gridLocal, &
    gridGlobal
@@ -150,7 +140,6 @@ subroutine Polarisation_init(temperature)
  integer(pInt) :: proc
  character(len=1024) :: rankStr
  
- call Utilities_init()
  if (worldrank == 0_pInt) then
    write(6,'(/,a)') ' <<<+-  DAMASK_spectral_solverPolarisation init  -+>>>'
    write(6,'(a)') ' $Id$'
@@ -169,6 +158,7 @@ subroutine Polarisation_init(temperature)
 !--------------------------------------------------------------------------------------------------
 ! PETSc Init
  call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr)
+ call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) 
  allocate(localK(worldsize), source = 0); localK(worldrank+1) = gridLocal(3)
  do proc = 1, worldsize
    call MPI_Bcast(localK(proc),1,MPI_INTEGER,proc-1,PETSC_COMM_WORLD,ierr)
@@ -182,10 +172,10 @@ subroutine Polarisation_init(temperature)
         gridLocal (1),gridLocal (2),localK, &                                                       ! local grid
         da,ierr)                                                                                    ! handle, error
  CHKERRQ(ierr)
+ call SNESSetDM(snes,da,ierr); CHKERRQ(ierr)
  call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr)
  call DMDASNESSetFunctionLocal(da,INSERT_VALUES,Polarisation_formResidual,dummy,ierr)
  CHKERRQ(ierr)
- call SNESSetDM(snes,da,ierr); CHKERRQ(ierr)
  call SNESSetConvergenceTest(snes,Polarisation_converged,dummy,PETSC_NULL_FUNCTION,ierr)
  CHKERRQ(ierr)
  call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr)
@@ -265,7 +255,7 @@ end subroutine Polarisation_init
 !--------------------------------------------------------------------------------------------------
 type(tSolutionState) function &
   Polarisation_solution(incInfoIn,guess,timeinc,timeinc_old,loadCaseTime,P_BC,F_BC,temperature_bc, &
-                                                                               rotation_BC,density)
+                                                                               rotation_BC)
  use numerics, only: &
    update_gamma
  use math, only: &
@@ -277,8 +267,6 @@ type(tSolutionState) function &
  use FEsolving, only: &
    restartWrite, &
    terminallyIll
- use numerics, only: &
-   worldrank
  
  implicit none
 
@@ -288,8 +276,7 @@ type(tSolutionState) function &
    timeinc, &                                                                                       !< increment in time for current solution
    timeinc_old, &                                                                                   !< increment in time of last increment
    loadCaseTime, &                                                                                  !< remaining time of current load case
-   temperature_bc, &
+   temperature_bc
-   density
  logical, intent(in) :: &
    guess
  type(tBoundaryCondition),      intent(in) :: &
@@ -324,7 +311,6 @@ type(tSolutionState) function &
  params%timeinc = timeinc
  params%timeincOld = timeinc_old
  params%temperature = temperature_bc
- params%density = density
 
 !--------------------------------------------------------------------------------------------------
 ! solve BVP 
@@ -363,16 +349,14 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
    math_transpose33, &
    math_mul3333xx33, &
    math_invSym3333, &
-   math_mul33x33, &
+   math_mul33x33
-   PI
  use DAMASK_spectral_Utilities, only: &
    wgt, &
-   field_realMPI, &
+   tensorField_realMPI, &
-   field_fourierMPI, &
+   utilities_FFTtensorForward, &
-   Utilities_FFTforward, &
+   utilities_fourierGammaConvolution, &
-   Utilities_fourierConvolution, &
    Utilities_inverseLaplace, &
-   Utilities_FFTbackward, &
+   utilities_FFTtensorBackward, &
    Utilities_constitutiveResponse, &
    Utilities_divergenceRMS, &
    Utilities_curlRMS
@@ -444,9 +428,9 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
 
 !--------------------------------------------------------------------------------------------------
 ! 
- field_realMPI = 0.0_pReal
+ tensorField_realMPI = 0.0_pReal
  do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt, gridLocal(1)
-   field_realMPI(1:3,1:3,i,j,k) = &
+   tensorField_realMPI(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))
@@ -454,13 +438,13 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
  
 !--------------------------------------------------------------------------------------------------
 ! doing convolution in Fourier space 
- call Utilities_FFTforward()
+ call utilities_FFTtensorForward()
- call Utilities_fourierConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC)) 
+ call utilities_fourierGammaConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC)) 
- call Utilities_FFTbackward()
+ call utilities_FFTtensorBackward()
 
 !--------------------------------------------------------------------------------------------------
 ! constructing residual                         
- residual_F_tau = polarBeta*F - field_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3))
+ residual_F_tau = polarBeta*F - tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3))
 
 !--------------------------------------------------------------------------------------------------
 ! evaluate constitutive response
@@ -472,11 +456,11 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
 
 !--------------------------------------------------------------------------------------------------
 ! calculate divergence
- field_realMPI = 0.0_pReal
+ tensorField_realMPI = 0.0_pReal
- field_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = residual_F
+ tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = residual_F
- call Utilities_FFTforward()
+ call utilities_FFTtensorForward()
  err_div = Utilities_divergenceRMS()
- call Utilities_FFTbackward()
+ call utilities_FFTtensorBackward()
  
 !--------------------------------------------------------------------------------------------------
 ! constructing residual
@@ -492,11 +476,11 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
  
 !--------------------------------------------------------------------------------------------------
 ! calculating curl
- field_realMPI = 0.0_pReal
+ tensorField_realMPI = 0.0_pReal
- field_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = F
+ tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = F
- call Utilities_FFTforward()
+ call utilities_FFTtensorForward()
  err_curl = Utilities_curlRMS()
- call Utilities_FFTbackward()
+ call utilities_FFTtensorBackward()
  
 end subroutine Polarisation_formResidual
 
@@ -727,7 +711,6 @@ subroutine Polarisation_destroy()
  call VecDestroy(solution_vec,ierr); CHKERRQ(ierr)
  call SNESDestroy(snes,ierr); CHKERRQ(ierr)
  call DMDestroy(da,ierr); CHKERRQ(ierr)
- call Utilities_destroy()
 
 end subroutine Polarisation_destroy
 

code/constitutive.f90

@@ -950,7 +950,7 @@ end subroutine constitutive_collectDotState
 !> @brief for constitutive models having an instantaneous change of state (so far, only nonlocal)
 !> will return false if delta state is not needed/supported by the constitutive model
 !--------------------------------------------------------------------------------------------------
-logical function constitutive_collectDeltaState(Tstar_v, Fe, ipc, ip, el)
+subroutine constitutive_collectDeltaState(Tstar_v, Fe, ipc, ip, el)
  use prec, only: &
    pReal, &
    pLongInt
@@ -998,26 +998,19 @@ logical function constitutive_collectDeltaState(Tstar_v, Fe, ipc, ip, el)
    call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
 
  select case (phase_plasticity(material_phase(ipc,ip,el)))
-
    case (PLASTICITY_NONLOCAL_ID)
-     constitutive_collectDeltaState = .true.
      call plastic_nonlocal_deltaState(Tstar_v,ip,el)
-   case default
-     constitutive_collectDeltaState = .false.
 
  end select
 
  do mySource = 1_pInt, phase_Nsources(material_phase(ipc,ip,el))
    select case (phase_source(mySource,material_phase(ipc,ip,el)))
      case (SOURCE_damage_isoBrittle_ID)
-       constitutive_collectDeltaState = constitutive_collectDeltaState .and. .true.
        call source_damage_isoBrittle_deltaState  (constitutive_homogenizedC(ipc,ip,el), Fe, &
                                                   ipc, ip, el)
      case (SOURCE_vacancy_irradiation_ID)
-       constitutive_collectDeltaState = constitutive_collectDeltaState .and. .true.
        call source_vacancy_irradiation_deltaState(ipc, ip, el)
      case (SOURCE_vacancy_thermalfluc_ID)
-       constitutive_collectDeltaState = constitutive_collectDeltaState .and. .true.
        call source_vacancy_thermalfluc_deltaState(ipc, ip, el)
 
    end select
@@ -1033,7 +1026,7 @@ logical function constitutive_collectDeltaState(Tstar_v, Fe, ipc, ip, el)
    !$OMP END CRITICAL (debugTimingDeltaState)
  endif
 
-end function constitutive_collectDeltaState
+end subroutine constitutive_collectDeltaState
 
 
 !--------------------------------------------------------------------------------------------------

code/crystallite.f90

@@ -3337,6 +3337,8 @@ logical function crystallite_stateJump(g,i,e)
    debug_g
  use material, only: &
    plasticState, &
+   sourceState, &
+   phase_Nsources, &
    mappingConstitutive
  use constitutive, only: &
    constitutive_collectDeltaState
@@ -3350,32 +3352,41 @@ logical function crystallite_stateJump(g,i,e)
  integer(pInt) :: &
    c, &
    p, &
-   mySizePlasticDotState
+   mySource, &
+   mySizePlasticDeltaState, &
+   mySizeSourceDeltaState
 
  c = mappingConstitutive(1,g,i,e)
  p = mappingConstitutive(2,g,i,e)
- if (constitutive_collectDeltaState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe(1:3,1:3,g,i,e), g,i,e)) then
+ call constitutive_collectDeltaState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe(1:3,1:3,g,i,e), g,i,e)
-   mySizePlasticDotState = plasticState(p)%sizeDotState
+ mySizePlasticDeltaState = plasticState(p)%sizeDeltaState
-   if( any(plasticState(p)%deltaState(:,c) /= plasticState(p)%deltaState(:,c))) then    ! NaN occured in deltaState
+ if( any(plasticState(p)%deltaState(:,c) /= plasticState(p)%deltaState(:,c))) then    ! NaN occured in deltaState
+   crystallite_stateJump = .false.
+   return
+ endif
+ plasticState(p)%state(1:mySizePlasticDeltaState,c) = plasticState(p)%state(1:mySizePlasticDeltaState,c) + &
+                                                      plasticState(p)%deltaState(1:mySizePlasticDeltaState,c)
+ do mySource = 1_pInt, phase_Nsources(p)
+   mySizeSourceDeltaState = sourceState(p)%p(mySource)%sizeDeltaState
+   if( any(sourceState(p)%p(mySource)%deltaState(:,c) /= sourceState(p)%p(mySource)%deltaState(:,c))) then    ! NaN occured in deltaState
      crystallite_stateJump = .false.
      return
    endif
-   plasticState(p)%state(1:mySizePlasticDotState,c) = plasticState(p)%state(1:mySizePlasticDotState,c) + &
+   sourceState(p)%p(mySource)%state(1:mySizeSourceDeltaState,c) = &
-                                                      plasticState(p)%deltaState(1:mySizePlasticDotState,c)
+     sourceState(p)%p(mySource)%state(1:mySizeSourceDeltaState,c) + &
+     sourceState(p)%p(mySource)%deltaState(1:mySizeSourceDeltaState,c)
+ enddo
 
 #ifndef _OPENMP
-   p = mappingConstitutive(2,g,i,e)
+ if (any(plasticState(p)%deltaState(1:mySizePlasticDeltaState,c) /= 0.0_pReal) &
-   c = mappingConstitutive(1,g,i,e)
+     .and. iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
-   if (any(plasticState(p)%deltaState(1:mySizePlasticDotState,c) /= 0.0_pReal) &
+     .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
-       .and. iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
+             .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
-       .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
+   write(6,'(a,i8,1x,i2,1x,i3, /)') '<< CRYST >> update state at el ip g ',e,i,g
-               .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
+   write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> deltaState', plasticState(p)%deltaState(1:mySizePlasticDeltaState,c)
-     write(6,'(a,i8,1x,i2,1x,i3, /)') '<< CRYST >> update state at el ip g ',e,i,g
+   write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> new state',  plasticState(p)%state     (1:mySizePlasticDeltaState,c)
-     write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> deltaState', plasticState(p)%deltaState(1:mySizePlasticDotState,c)
-     write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> new state',  plasticState(p)%state     (1:mySizePlasticDotState,c)
-   endif
-#endif
  endif
+#endif
 
  crystallite_stateJump = .true.
 

code/damage_nonlocal.f90

@@ -253,6 +253,8 @@ end subroutine damage_nonlocal_getSourceAndItsTangent
 !> @brief returns homogenized non local damage diffusion tensor in reference configuration
 !--------------------------------------------------------------------------------------------------
 function damage_nonlocal_getDiffusion33(ip,el)
+ use numerics, only: &
+   charLength
  use lattice, only: &
    lattice_DamageDiffusion33
  use material, only: &
@@ -280,6 +282,7 @@ function damage_nonlocal_getDiffusion33(ip,el)
  enddo
 
  damage_nonlocal_getDiffusion33 = &
+   charLength*charLength* &
    damage_nonlocal_getDiffusion33/ &
    homogenization_Ngrains(homog)
  

code/plastic_disloKMC.f90

@@ -201,7 +201,7 @@ subroutine plastic_disloKMC_init(fileUnit)
                   Nchunks_TwinSlip = 0_pInt, Nchunks_TwinTwin = 0_pInt, &
                   Nchunks_SlipFamilies = 0_pInt, Nchunks_TwinFamilies = 0_pInt, Nchunks_nonSchmid = 0_pInt, &
                   offset_slip, index_myFamily, index_otherFamily
- integer(pInt) :: sizeState, sizeDotState
+ integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
  integer(pInt) :: NofMyPhase
  character(len=65536) :: &
    tag  = '', &
@@ -658,12 +658,14 @@ subroutine plastic_disloKMC_init(fileUnit)
 ! allocate state arrays
      sizeDotState              =   int(size(plastic_disloKMC_listBasicSlipStates),pInt) * ns &
                                  + int(size(plastic_disloKMC_listBasicTwinStates),pInt) * nt
+     sizeDeltaState            =   sizeDotState
      sizeState                 =   sizeDotState &
                                  + int(size(plastic_disloKMC_listDependentSlipStates),pInt) * ns &
                                  + int(size(plastic_disloKMC_listDependentTwinStates),pInt) * nt
                 
      plasticState(phase)%sizeState = sizeState
      plasticState(phase)%sizeDotState = sizeDotState
+     plasticState(phase)%sizeDeltaState = sizeDeltaState
      plasticState(phase)%sizePostResults = plastic_disloKMC_sizePostResults(instance)
      plasticState(phase)%nSlip = plastic_disloKMC_totalNslip(instance)
      plasticState(phase)%nTwin = 0_pInt
@@ -676,7 +678,7 @@ subroutine plastic_disloKMC_init(fileUnit)
      allocate(plasticState(phase)%state_backup        (sizeState,NofMyPhase),     source=0.0_pReal)
 
      allocate(plasticState(phase)%dotState            (sizeDotState,NofMyPhase),  source=0.0_pReal)
-     allocate(plasticState(phase)%deltaState          (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     allocate(plasticState(phase)%deltaState        (sizeDeltaState,NofMyPhase),  source=0.0_pReal)
      allocate(plasticState(phase)%dotState_backup     (sizeDotState,NofMyPhase),  source=0.0_pReal)
      if (any(numerics_integrator == 1_pInt)) then
        allocate(plasticState(phase)%previousDotState  (sizeDotState,NofMyPhase),  source=0.0_pReal)

code/plastic_disloUCLA.f90

@@ -206,7 +206,7 @@ subroutine plastic_disloUCLA_init(fileUnit)
                   Nchunks_TwinSlip = 0_pInt, Nchunks_TwinTwin = 0_pInt, &
                   Nchunks_SlipFamilies = 0_pInt, Nchunks_TwinFamilies = 0_pInt, Nchunks_nonSchmid = 0_pInt, &
                   offset_slip, index_myFamily, index_otherFamily
- integer(pInt) :: sizeState, sizeDotState
+ integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
  integer(pInt) :: NofMyPhase
  character(len=65536) :: &
    tag  = '', &
@@ -677,12 +677,14 @@ subroutine plastic_disloUCLA_init(fileUnit)
 ! allocate state arrays
      sizeDotState              =   int(size(plastic_disloUCLA_listBasicSlipStates),pInt) * ns &
                                  + int(size(plastic_disloUCLA_listBasicTwinStates),pInt) * nt
+     sizeDeltaState            =   0_pInt
      sizeState                 =   sizeDotState &
                                  + int(size(plastic_disloUCLA_listDependentSlipStates),pInt) * ns &
                                  + int(size(plastic_disloUCLA_listDependentTwinStates),pInt) * nt
                 
      plasticState(phase)%sizeState = sizeState
      plasticState(phase)%sizeDotState = sizeDotState
+     plasticState(phase)%sizeDeltaState = sizeDeltaState
      plasticState(phase)%sizePostResults = plastic_disloUCLA_sizePostResults(instance)
      plasticState(phase)%nSlip = plastic_disloucla_totalNslip(instance)
      plasticState(phase)%nTwin = plastic_disloucla_totalNtwin(instance)
@@ -695,7 +697,7 @@ subroutine plastic_disloUCLA_init(fileUnit)
      allocate(plasticState(phase)%state_backup        (sizeState,NofMyPhase),     source=0.0_pReal)
 
      allocate(plasticState(phase)%dotState            (sizeDotState,NofMyPhase),  source=0.0_pReal)
-     allocate(plasticState(phase)%deltaState          (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     allocate(plasticState(phase)%deltaState        (sizeDeltaState,NofMyPhase),  source=0.0_pReal)
      allocate(plasticState(phase)%dotState_backup     (sizeDotState,NofMyPhase),  source=0.0_pReal)
      if (any(numerics_integrator == 1_pInt)) then
        allocate(plasticState(phase)%previousDotState  (sizeDotState,NofMyPhase),  source=0.0_pReal)

code/plastic_dislotwin.f90

@@ -233,7 +233,7 @@ subroutine plastic_dislotwin_init(fileUnit)
                   Nchunks_TwinSlip = 0_pInt, Nchunks_TwinTwin = 0_pInt, &
                   Nchunks_SlipFamilies = 0_pInt, Nchunks_TwinFamilies = 0_pInt, Nchunks_TransFamilies = 0_pInt, &
                   offset_slip, index_myFamily, index_otherFamily
- integer(pInt) :: sizeState, sizeDotState
+ integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
  integer(pInt) :: NofMyPhase   
  character(len=65536) :: &
    tag  = '', &
@@ -822,12 +822,14 @@ subroutine plastic_dislotwin_init(fileUnit)
      sizeDotState              =   int(size(plastic_dislotwin_listBasicSlipStates),pInt) * ns &
                                  + int(size(plastic_dislotwin_listBasicTwinStates),pInt) * nt &
                                  + int(size(plastic_dislotwin_listBasicTransStates),pInt) * nr
+     sizeDeltaState            =   0_pInt
      sizeState                 =   sizeDotState &
                                  + int(size(plastic_dislotwin_listDependentSlipStates),pInt) * ns &
                                  + int(size(plastic_dislotwin_listDependentTwinStates),pInt) * nt
                 
      plasticState(phase)%sizeState = sizeState
      plasticState(phase)%sizeDotState = sizeDotState
+     plasticState(phase)%sizeDeltaState = sizeDeltaState
      plasticState(phase)%sizePostResults = plastic_dislotwin_sizePostResults(instance)
      plasticState(phase)%nSlip = plastic_dislotwin_totalNslip(instance)
      plasticState(phase)%nTwin = plastic_dislotwin_totalNtwin(instance)
@@ -840,7 +842,7 @@ subroutine plastic_dislotwin_init(fileUnit)
      allocate(plasticState(phase)%state_backup        (sizeState,NofMyPhase),     source=0.0_pReal)
 
      allocate(plasticState(phase)%dotState            (sizeDotState,NofMyPhase),  source=0.0_pReal)
-     allocate(plasticState(phase)%deltaState          (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     allocate(plasticState(phase)%deltaState        (sizeDeltaState,NofMyPhase),  source=0.0_pReal)
      allocate(plasticState(phase)%dotState_backup     (sizeDotState,NofMyPhase),  source=0.0_pReal)
      if (any(numerics_integrator == 1_pInt)) then
        allocate(plasticState(phase)%previousDotState  (sizeDotState,NofMyPhase),  source=0.0_pReal)

code/plastic_j2.f90

@@ -141,7 +141,8 @@ subroutine plastic_j2_init(fileUnit)
    instance, &
    mySize, &
    sizeDotState, &
-   sizeState
+   sizeState, &
+   sizeDeltaState
  character(len=65536) :: &
    tag  = '', &
    line = ''
@@ -324,8 +325,10 @@ subroutine plastic_j2_init(fileUnit)
 ! allocate state arrays
      sizeState    = 2_pInt
      sizeDotState = sizeState
+     sizeDeltaState = 0_pInt
      plasticState(phase)%sizeState = sizeState
      plasticState(phase)%sizeDotState = sizeDotState
+     plasticState(phase)%sizeDeltaState = sizeDeltaState
      plasticState(phase)%sizePostResults = plastic_j2_sizePostResults(instance)
      plasticState(phase)%nSlip = 1
      plasticState(phase)%nTwin = 0
@@ -341,6 +344,7 @@ subroutine plastic_j2_init(fileUnit)
      allocate(plasticState(phase)%state              (   sizeState,NofMyPhase),source=0.0_pReal)
      allocate(plasticState(phase)%state_backup       (   sizeState,NofMyPhase),source=0.0_pReal)
      allocate(plasticState(phase)%dotState           (sizeDotState,NofMyPhase),source=0.0_pReal)
+     allocate(plasticState(phase)%deltaState        (sizeDeltaState,NofMyPhase),  source=0.0_pReal)
      allocate(plasticState(phase)%dotState_backup    (sizeDotState,NofMyPhase),source=0.0_pReal)
      if (any(numerics_integrator == 1_pInt)) then
        allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase),source=0.0_pReal)

code/plastic_none.f90

@@ -52,7 +52,8 @@ subroutine plastic_none_init
    phase, &
    NofMyPhase, &
    sizeState, &
-   sizeDotState
+   sizeDotState, &
+   sizeDeltaState
  
  mainProcess: if (worldrank == 0) then 
    write(6,'(/,a)')   ' <<<+-  constitutive_'//PLASTICITY_NONE_label//' init  -+>>>'
@@ -75,6 +76,8 @@ subroutine plastic_none_init
      plasticState(phase)%sizeState = sizeState
      sizeDotState = sizeState
      plasticState(phase)%sizeDotState = sizeDotState
+     sizeDeltaState = 0_pInt
+     plasticState(phase)%sizeDeltaState = sizeDeltaState
      plasticState(phase)%sizePostResults = 0_pInt
      plasticState(phase)%nSlip  = 0_pInt
      plasticState(phase)%nTwin  = 0_pInt
@@ -87,6 +90,7 @@ subroutine plastic_none_init
      allocate(plasticState(phase)%state_backup       (sizeState,NofMyPhase))
 
      allocate(plasticState(phase)%dotState           (sizeDotState,NofMyPhase))
+     allocate(plasticState(phase)%deltaState        (sizeDeltaState,NofMyPhase))
      allocate(plasticState(phase)%dotState_backup    (sizeDotState,NofMyPhase))
      if (any(numerics_integrator == 1_pInt)) then
        allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase))

code/plastic_nonlocal.f90

@@ -330,7 +330,7 @@ integer(pInt)          ::                   phase, &
    tag  = '', &
    line = ''
 
- integer(pInt) :: sizeState, sizeDotState,sizeDependentState
+ integer(pInt) :: sizeState, sizeDotState,sizeDependentState, sizeDeltaState
 
 
  integer(pInt) :: NofMyPhase 
@@ -1144,6 +1144,7 @@ allocate(nonSchmidProjection(3,3,4,maxTotalNslip,maxNinstances),
      sizeDependentState        = int(size(DEPENDENTSTATES),pInt) * ns
      sizeState                 = sizeDotState + sizeDependentState &
                                + int(size(OTHERSTATES),pInt) * ns
+     sizeDeltaState            = sizeDotState
 
      !*** determine indices to state array
 
@@ -1299,6 +1300,7 @@ allocate(nonSchmidProjection(3,3,4,maxTotalNslip,maxNinstances),
                 
      plasticState(phase)%sizeState    = sizeState
      plasticState(phase)%sizeDotState = sizeDotState
+     plasticState(phase)%sizeDeltaState = sizeDeltaState
      plasticState(phase)%sizePostResults = plastic_nonlocal_sizePostResults(instance)
      plasticState(phase)%nonlocal = .true.
      plasticState(phase)%nSlip = totalNslip(instance)
@@ -1312,7 +1314,7 @@ allocate(nonSchmidProjection(3,3,4,maxTotalNslip,maxNinstances),
      allocate(plasticState(phase)%state_backup        (sizeState,NofMyPhase),     source=0.0_pReal)
 
      allocate(plasticState(phase)%dotState            (sizeDotState,NofMyPhase),  source=0.0_pReal)
-     allocate(plasticState(phase)%deltaState          (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     allocate(plasticState(phase)%deltaState        (sizeDeltaState,NofMyPhase),  source=0.0_pReal)
      allocate(plasticState(phase)%dotState_backup     (sizeDotState,NofMyPhase),  source=0.0_pReal)
      if (any(numerics_integrator == 1_pInt)) then
        allocate(plasticState(phase)%previousDotState  (sizeDotState,NofMyPhase),  source=0.0_pReal)

code/plastic_phenopowerlaw.f90

@@ -158,7 +158,7 @@ subroutine plastic_phenopowerlaw_init(fileUnit)
    Nchunks_TransFamilies = 0_pInt, Nchunks_nonSchmid = 0_pInt, &
    NipcMyPhase, &
    offset_slip, index_myFamily, index_otherFamily, &
-   mySize=0_pInt,sizeState,sizeDotState
+   mySize=0_pInt,sizeState,sizeDotState, sizeDeltaState
  character(len=65536) :: &
    tag  = '', &
    line = ''
@@ -552,8 +552,10 @@ subroutine plastic_phenopowerlaw_init(fileUnit)
                + plastic_phenopowerlaw_totalNtwin(instance)                           ! accshear_twin
                
      sizeDotState = sizeState
+     sizeDeltaState = 0_pInt
      plasticState(phase)%sizeState = sizeState
      plasticState(phase)%sizeDotState = sizeDotState
+     plasticState(phase)%sizeDeltaState = sizeDeltaState
      plasticState(phase)%sizePostResults = plastic_phenopowerlaw_sizePostResults(instance)
      plasticState(phase)%nSlip =plastic_phenopowerlaw_totalNslip(instance)
      plasticState(phase)%nTwin =plastic_phenopowerlaw_totalNtwin(instance)
@@ -565,6 +567,7 @@ subroutine plastic_phenopowerlaw_init(fileUnit)
      allocate(plasticState(phase)%state              (   sizeState,NipcMyPhase), source=0.0_pReal)
      allocate(plasticState(phase)%state_backup       (   sizeState,NipcMyPhase), source=0.0_pReal)
      allocate(plasticState(phase)%dotState           (sizeDotState,NipcMyPhase), source=0.0_pReal)
+     allocate(plasticState(phase)%deltaState       (sizeDeltaState,NipcMyPhase), source=0.0_pReal)
      allocate(plasticState(phase)%dotState_backup    (sizeDotState,NipcMyPhase), source=0.0_pReal)
      if (any(numerics_integrator == 1_pInt)) then
        allocate(plasticState(phase)%previousDotState (sizeDotState,NipcMyPhase),source=0.0_pReal)

code/plastic_titanmod.f90

@@ -237,7 +237,7 @@ subroutine plastic_titanmod_init(fileUnit)
    Nchunks_SlipFamilies = 0_pInt, Nchunks_TwinFamilies = 0_pInt, &
    offset_slip, mySize, &
    maxTotalNslip,maxTotalNtwin, maxNinstance
- integer(pInt) :: sizeState, sizeDotState
+ integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
  integer(pInt) :: NofMyPhase = 0_pInt  
  character(len=65536) :: &
    tag  = '', &
@@ -812,6 +812,7 @@ subroutine plastic_titanmod_init(fileUnit)
      sizeState    = sizeDotState+ &
        size(plastic_titanmod_listDependentSlipStates)*ns + &
        size(plastic_titanmod_listDependentTwinStates)*nt
+     sizeDeltaState = 0_pInt
 
 !--------------------------------------------------------------------------------------------------
 ! determine size of postResults array  
@@ -848,6 +849,7 @@ subroutine plastic_titanmod_init(fileUnit)
 ! Determine size of state array                
      plasticState(phase)%sizeState = sizeState
      plasticState(phase)%sizeDotState = sizeDotState
+     plasticState(phase)%sizeDeltaState = sizeDeltaState
      plasticState(phase)%sizePostResults = plastic_titanmod_sizePostResults(instance)
      plasticState(phase)%nSlip =plastic_titanmod_totalNslip(instance)
      plasticState(phase)%nTwin = 0_pInt
@@ -860,7 +862,7 @@ subroutine plastic_titanmod_init(fileUnit)
      allocate(plasticState(phase)%state_backup        (sizeState,NofMyPhase),     source=0.0_pReal)
 
      allocate(plasticState(phase)%dotState            (sizeDotState,NofMyPhase),  source=0.0_pReal)
-     allocate(plasticState(phase)%deltaState          (sizeDotState,NofMyPhase),     source=0.0_pReal)
+     allocate(plasticState(phase)%deltaState        (sizeDeltaState,NofMyPhase),  source=0.0_pReal)
      allocate(plasticState(phase)%dotState_backup     (sizeDotState,NofMyPhase),  source=0.0_pReal)
      if (any(numerics_integrator == 1_pInt)) then
        allocate(plasticState(phase)%previousDotState  (sizeDotState,NofMyPhase),  source=0.0_pReal)

code/prec.f90

@@ -70,6 +70,7 @@ use ifport
    integer(pInt) :: &
      sizeState = 0_pInt , &                                                                         !< size of state
      sizeDotState = 0_pInt, &                                                                       !< size of dot state, i.e. parts of the state that are integrated
+     sizeDeltaState = 0_pInt, &                                                                     !< size of delta state, i.e. parts of the state that have discontinuous rates
      sizePostResults = 0_pInt                                                                       !< size of output data
    real(pReal), allocatable, dimension(:) :: &
      atolState

code/source_damage_anisoBrittle.f90

@@ -106,7 +106,7 @@ subroutine source_damage_anisoBrittle_init(fileUnit)
  integer(pInt), parameter :: MAXNCHUNKS = 7_pInt
  integer(pInt), dimension(1+2*MAXNCHUNKS) :: positions
  integer(pInt) :: maxNinstance,mySize=0_pInt,phase,instance,source,sourceOffset,o
- integer(pInt) :: sizeState, sizeDotState
+ integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
  integer(pInt) :: NofMyPhase   
  integer(pInt) :: Nchunks_CleavageFamilies = 0_pInt, j   
  character(len=65536) :: &
@@ -255,10 +255,12 @@ subroutine source_damage_anisoBrittle_init(fileUnit)
 !--------------------------------------------------------------------------------------------------
 ! Determine size of state array
      sizeDotState              = 1_pInt
+     sizeDeltaState            = 0_pInt
      sizeState                 = 1_pInt
 
-     sourceState(phase)%p(sourceOffset)%sizeState = sizeState
+     sourceState(phase)%p(sourceOffset)%sizeState      = sizeState
-     sourceState(phase)%p(sourceOffset)%sizeDotState = sizeDotState
+     sourceState(phase)%p(sourceOffset)%sizeDotState   = sizeDotState
+     sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState
      sourceState(phase)%p(sourceOffset)%sizePostResults = source_damage_anisoBrittle_sizePostResults(instance)
      allocate(sourceState(phase)%p(sourceOffset)%aTolState           (sizeState),                &
               source=source_damage_anisoBrittle_aTol(instance))
@@ -269,7 +271,7 @@ subroutine source_damage_anisoBrittle_init(fileUnit)
      allocate(sourceState(phase)%p(sourceOffset)%state_backup        (sizeState,NofMyPhase),     source=0.0_pReal)
 
      allocate(sourceState(phase)%p(sourceOffset)%dotState            (sizeDotState,NofMyPhase),  source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%deltaState          (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%deltaState        (sizeDeltaState,NofMyPhase),  source=0.0_pReal)
      allocate(sourceState(phase)%p(sourceOffset)%dotState_backup     (sizeDotState,NofMyPhase),  source=0.0_pReal)
      if (any(numerics_integrator == 1_pInt)) then
        allocate(sourceState(phase)%p(sourceOffset)%previousDotState  (sizeDotState,NofMyPhase),  source=0.0_pReal)

code/source_damage_anisoDuctile.f90

@@ -110,7 +110,7 @@ subroutine source_damage_anisoDuctile_init(fileUnit)
  integer(pInt), parameter :: MAXNCHUNKS = 7_pInt
  integer(pInt), dimension(1+2*MAXNCHUNKS) :: positions
  integer(pInt) :: maxNinstance,mySize=0_pInt,phase,instance,source,sourceOffset,o
- integer(pInt) :: sizeState, sizeDotState
+ integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
  integer(pInt) :: NofMyPhase   
  integer(pInt) :: Nchunks_SlipFamilies = 0_pInt, j   
  character(len=65536) :: &
@@ -258,9 +258,11 @@ subroutine source_damage_anisoDuctile_init(fileUnit)
 !--------------------------------------------------------------------------------------------------
 ! Determine size of state array
      sizeDotState              = 1_pInt
+     sizeDeltaState            = 0_pInt
      sizeState                 = 1_pInt
-     sourceState(phase)%p(sourceOffset)%sizeState = sizeState
+     sourceState(phase)%p(sourceOffset)%sizeState      = sizeState
-     sourceState(phase)%p(sourceOffset)%sizeDotState = sizeDotState
+     sourceState(phase)%p(sourceOffset)%sizeDotState   = sizeDotState
+     sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState
      sourceState(phase)%p(sourceOffset)%sizePostResults = source_damage_anisoDuctile_sizePostResults(instance)
      allocate(sourceState(phase)%p(sourceOffset)%aTolState           (sizeState),                &
               source=source_damage_anisoDuctile_aTol(instance))
@@ -271,7 +273,7 @@ subroutine source_damage_anisoDuctile_init(fileUnit)
      allocate(sourceState(phase)%p(sourceOffset)%state_backup        (sizeState,NofMyPhase),     source=0.0_pReal)
 
      allocate(sourceState(phase)%p(sourceOffset)%dotState            (sizeDotState,NofMyPhase),  source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%deltaState          (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%deltaState        (sizeDeltaState,NofMyPhase),  source=0.0_pReal)
      allocate(sourceState(phase)%p(sourceOffset)%dotState_backup     (sizeDotState,NofMyPhase),  source=0.0_pReal)
      if (any(numerics_integrator == 1_pInt)) then
        allocate(sourceState(phase)%p(sourceOffset)%previousDotState  (sizeDotState,NofMyPhase),  source=0.0_pReal)

code/source_damage_isoBrittle.f90

@@ -92,7 +92,7 @@ subroutine source_damage_isoBrittle_init(fileUnit)
  integer(pInt), parameter :: MAXNCHUNKS = 7_pInt
  integer(pInt), dimension(1+2*MAXNCHUNKS) :: positions
  integer(pInt) :: maxNinstance,mySize=0_pInt,phase,instance,source,sourceOffset,o
- integer(pInt) :: sizeState, sizeDotState
+ integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
  integer(pInt) :: NofMyPhase   
  character(len=65536) :: &
    tag  = '', &
@@ -204,10 +204,12 @@ subroutine source_damage_isoBrittle_init(fileUnit)
      enddo outputsLoop
 ! Determine size of state array
      sizeDotState              =   1_pInt
+     sizeDeltaState            =   1_pInt
      sizeState                 =   1_pInt
                 
-     sourceState(phase)%p(sourceOffset)%sizeState = sizeState
+     sourceState(phase)%p(sourceOffset)%sizeState      = sizeState
-     sourceState(phase)%p(sourceOffset)%sizeDotState = sizeDotState
+     sourceState(phase)%p(sourceOffset)%sizeDotState   = sizeDotState
+     sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState
      sourceState(phase)%p(sourceOffset)%sizePostResults = source_damage_isoBrittle_sizePostResults(instance)
      allocate(sourceState(phase)%p(sourceOffset)%aTolState           (sizeState),                &
               source=source_damage_isoBrittle_aTol(instance))
@@ -218,7 +220,7 @@ subroutine source_damage_isoBrittle_init(fileUnit)
      allocate(sourceState(phase)%p(sourceOffset)%state_backup        (sizeState,NofMyPhase),     source=0.0_pReal)
 
      allocate(sourceState(phase)%p(sourceOffset)%dotState            (sizeDotState,NofMyPhase),  source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%deltaState          (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%deltaState        (sizeDeltaState,NofMyPhase),  source=0.0_pReal)
      allocate(sourceState(phase)%p(sourceOffset)%dotState_backup     (sizeDotState,NofMyPhase),  source=0.0_pReal)
      if (any(numerics_integrator == 1_pInt)) then
        allocate(sourceState(phase)%p(sourceOffset)%previousDotState  (sizeDotState,NofMyPhase),  source=0.0_pReal)

code/source_damage_isoDuctile.f90

@@ -93,7 +93,7 @@ subroutine source_damage_isoDuctile_init(fileUnit)
  integer(pInt), parameter :: MAXNCHUNKS = 7_pInt
  integer(pInt), dimension(1+2*MAXNCHUNKS) :: positions
  integer(pInt) :: maxNinstance,mySize=0_pInt,phase,instance,source,sourceOffset,o
- integer(pInt) :: sizeState, sizeDotState
+ integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
  integer(pInt) :: NofMyPhase   
  character(len=65536) :: &
    tag  = '', &
@@ -208,11 +208,13 @@ subroutine source_damage_isoDuctile_init(fileUnit)
        endif
      enddo outputsLoop
 ! Determine size of state array
-     sizeDotState              =   0_pInt
+     sizeDotState              =   1_pInt
+     sizeDeltaState            =   0_pInt
      sizeState                 =   1_pInt
                 
-     sourceState(phase)%p(sourceOffset)%sizeState = sizeState
+     sourceState(phase)%p(sourceOffset)%sizeState      = sizeState
-     sourceState(phase)%p(sourceOffset)%sizeDotState = sizeDotState
+     sourceState(phase)%p(sourceOffset)%sizeDotState   = sizeDotState
+     sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState
      sourceState(phase)%p(sourceOffset)%sizePostResults = source_damage_isoDuctile_sizePostResults(instance)
      allocate(sourceState(phase)%p(sourceOffset)%aTolState           (sizeState),                &
               source=source_damage_isoDuctile_aTol(instance))
@@ -223,7 +225,7 @@ subroutine source_damage_isoDuctile_init(fileUnit)
      allocate(sourceState(phase)%p(sourceOffset)%state_backup        (sizeState,NofMyPhase),     source=0.0_pReal)
 
      allocate(sourceState(phase)%p(sourceOffset)%dotState            (sizeDotState,NofMyPhase),  source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%deltaState          (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%deltaState        (sizeDeltaState,NofMyPhase),  source=0.0_pReal)
      allocate(sourceState(phase)%p(sourceOffset)%dotState_backup     (sizeDotState,NofMyPhase),  source=0.0_pReal)
      if (any(numerics_integrator == 1_pInt)) then
        allocate(sourceState(phase)%p(sourceOffset)%previousDotState  (sizeDotState,NofMyPhase),  source=0.0_pReal)

code/source_thermal_dissipation.f90

@@ -79,7 +79,7 @@ subroutine source_thermal_dissipation_init(fileUnit)
  integer(pInt), parameter :: MAXNCHUNKS = 7_pInt
  integer(pInt), dimension(1+2*MAXNCHUNKS) :: positions
  integer(pInt) :: maxNinstance,phase,instance,source,sourceOffset
- integer(pInt) :: sizeState, sizeDotState
+ integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
  integer(pInt) :: NofMyPhase   
  character(len=65536) :: &
    tag  = '', &
@@ -152,9 +152,11 @@ subroutine source_thermal_dissipation_init(fileUnit)
      sourceOffset = source_thermal_dissipation_offset(phase)
 
      sizeDotState              =   0_pInt
+     sizeDeltaState            =   0_pInt
      sizeState                 =   0_pInt
      sourceState(phase)%p(sourceOffset)%sizeState =       sizeState
      sourceState(phase)%p(sourceOffset)%sizeDotState =    sizeDotState
+     sourceState(phase)%p(sourceOffset)%sizeDeltaState =  sizeDeltaState
      sourceState(phase)%p(sourceOffset)%sizePostResults = source_thermal_dissipation_sizePostResults(instance)
      allocate(sourceState(phase)%p(sourceOffset)%aTolState           (sizeState),                source=0.0_pReal)
      allocate(sourceState(phase)%p(sourceOffset)%state0              (sizeState,NofMyPhase),     source=0.0_pReal)
@@ -164,7 +166,7 @@ subroutine source_thermal_dissipation_init(fileUnit)
      allocate(sourceState(phase)%p(sourceOffset)%state_backup        (sizeState,NofMyPhase),     source=0.0_pReal)
 
      allocate(sourceState(phase)%p(sourceOffset)%dotState            (sizeDotState,NofMyPhase),  source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%deltaState          (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%deltaState        (sizeDeltaState,NofMyPhase),  source=0.0_pReal)
      allocate(sourceState(phase)%p(sourceOffset)%dotState_backup     (sizeDotState,NofMyPhase),  source=0.0_pReal)
      if (any(numerics_integrator == 1_pInt)) then
        allocate(sourceState(phase)%p(sourceOffset)%previousDotState  (sizeDotState,NofMyPhase),  source=0.0_pReal)

code/source_vacancy_irradiation.f90

@@ -81,7 +81,7 @@ subroutine source_vacancy_irradiation_init(fileUnit)
  integer(pInt), parameter :: MAXNCHUNKS = 7_pInt
  integer(pInt), dimension(1+2*MAXNCHUNKS) :: positions
  integer(pInt) :: maxNinstance,phase,instance,source,sourceOffset
- integer(pInt) :: sizeState, sizeDotState
+ integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
  integer(pInt) :: NofMyPhase   
  character(len=65536) :: &
    tag  = '', &
@@ -158,9 +158,11 @@ subroutine source_vacancy_irradiation_init(fileUnit)
      sourceOffset = source_vacancy_irradiation_offset(phase)
 
      sizeDotState              =   2_pInt
+     sizeDeltaState            =   2_pInt
      sizeState                 =   2_pInt
      sourceState(phase)%p(sourceOffset)%sizeState =       sizeState
      sourceState(phase)%p(sourceOffset)%sizeDotState =    sizeDotState
+     sourceState(phase)%p(sourceOffset)%sizeDeltaState =  sizeDeltaState
      sourceState(phase)%p(sourceOffset)%sizePostResults = source_vacancy_irradiation_sizePostResults(instance)
      allocate(sourceState(phase)%p(sourceOffset)%aTolState           (sizeState),                source=0.1_pReal)
      allocate(sourceState(phase)%p(sourceOffset)%state0              (sizeState,NofMyPhase),     source=0.0_pReal)
@@ -170,7 +172,7 @@ subroutine source_vacancy_irradiation_init(fileUnit)
      allocate(sourceState(phase)%p(sourceOffset)%state_backup        (sizeState,NofMyPhase),     source=0.0_pReal)
 
      allocate(sourceState(phase)%p(sourceOffset)%dotState            (sizeDotState,NofMyPhase),  source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%deltaState          (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%deltaState        (sizeDeltaState,NofMyPhase),  source=0.0_pReal)
      allocate(sourceState(phase)%p(sourceOffset)%dotState_backup     (sizeDotState,NofMyPhase),  source=0.0_pReal)
      if (any(numerics_integrator == 1_pInt)) then
        allocate(sourceState(phase)%p(sourceOffset)%previousDotState  (sizeDotState,NofMyPhase),  source=0.0_pReal)

code/source_vacancy_phenoplasticity.f90

@@ -79,7 +79,7 @@ subroutine source_vacancy_phenoplasticity_init(fileUnit)
  integer(pInt), parameter :: MAXNCHUNKS = 7_pInt
  integer(pInt), dimension(1+2*MAXNCHUNKS) :: positions
  integer(pInt) :: maxNinstance,phase,instance,source,sourceOffset
- integer(pInt) :: sizeState, sizeDotState
+ integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
  integer(pInt) :: NofMyPhase   
  character(len=65536) :: &
    tag  = '', &
@@ -152,9 +152,11 @@ subroutine source_vacancy_phenoplasticity_init(fileUnit)
      sourceOffset = source_vacancy_phenoplasticity_offset(phase)
 
      sizeDotState              =   0_pInt
+     sizeDeltaState            =   0_pInt
      sizeState                 =   0_pInt
      sourceState(phase)%p(sourceOffset)%sizeState =       sizeState
      sourceState(phase)%p(sourceOffset)%sizeDotState =    sizeDotState
+     sourceState(phase)%p(sourceOffset)%sizeDeltaState =  sizeDeltaState
      sourceState(phase)%p(sourceOffset)%sizePostResults = source_vacancy_phenoplasticity_sizePostResults(instance)
      allocate(sourceState(phase)%p(sourceOffset)%aTolState           (sizeState),                source=0.0_pReal)
      allocate(sourceState(phase)%p(sourceOffset)%state0              (sizeState,NofMyPhase),     source=0.0_pReal)
@@ -164,7 +166,7 @@ subroutine source_vacancy_phenoplasticity_init(fileUnit)
      allocate(sourceState(phase)%p(sourceOffset)%state_backup        (sizeState,NofMyPhase),     source=0.0_pReal)
 
      allocate(sourceState(phase)%p(sourceOffset)%dotState            (sizeDotState,NofMyPhase),  source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%deltaState          (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%deltaState        (sizeDeltaState,NofMyPhase),  source=0.0_pReal)
      allocate(sourceState(phase)%p(sourceOffset)%dotState_backup     (sizeDotState,NofMyPhase),  source=0.0_pReal)
      if (any(numerics_integrator == 1_pInt)) then
        allocate(sourceState(phase)%p(sourceOffset)%previousDotState  (sizeDotState,NofMyPhase),  source=0.0_pReal)

code/source_vacancy_thermalfluc.f90

@@ -80,7 +80,7 @@ subroutine source_vacancy_thermalfluc_init(fileUnit)
  integer(pInt), parameter :: MAXNCHUNKS = 7_pInt
  integer(pInt), dimension(1+2*MAXNCHUNKS) :: positions
  integer(pInt) :: maxNinstance,phase,instance,source,sourceOffset
- integer(pInt) :: sizeState, sizeDotState
+ integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
  integer(pInt) :: NofMyPhase   
  character(len=65536) :: &
    tag  = '', &
@@ -153,9 +153,11 @@ subroutine source_vacancy_thermalfluc_init(fileUnit)
      sourceOffset = source_vacancy_thermalfluc_offset(phase)
 
      sizeDotState              =   1_pInt
+     sizeDeltaState            =   1_pInt
      sizeState                 =   1_pInt
      sourceState(phase)%p(sourceOffset)%sizeState =       sizeState
      sourceState(phase)%p(sourceOffset)%sizeDotState =    sizeDotState
+     sourceState(phase)%p(sourceOffset)%sizeDeltaState =  sizeDeltaState
      sourceState(phase)%p(sourceOffset)%sizePostResults = source_vacancy_thermalfluc_sizePostResults(instance)
      allocate(sourceState(phase)%p(sourceOffset)%aTolState           (sizeState),                source=0.0_pReal)
      allocate(sourceState(phase)%p(sourceOffset)%state0              (sizeState,NofMyPhase),     source=0.0_pReal)
@@ -165,7 +167,7 @@ subroutine source_vacancy_thermalfluc_init(fileUnit)
      allocate(sourceState(phase)%p(sourceOffset)%state_backup        (sizeState,NofMyPhase),     source=0.0_pReal)
 
      allocate(sourceState(phase)%p(sourceOffset)%dotState            (sizeDotState,NofMyPhase),  source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%deltaState          (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%deltaState        (sizeDeltaState,NofMyPhase),  source=0.0_pReal)
      allocate(sourceState(phase)%p(sourceOffset)%dotState_backup     (sizeDotState,NofMyPhase),  source=0.0_pReal)
      if (any(numerics_integrator == 1_pInt)) then
        allocate(sourceState(phase)%p(sourceOffset)%previousDotState  (sizeDotState,NofMyPhase),  source=0.0_pReal)