diff --git a/code/constitutive.f90 b/code/constitutive.f90 index 37804f5d1..ed4ea06e2 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -138,9 +138,7 @@ subroutine constitutive_init use constitutive_phenopowerlaw use constitutive_dislotwin use constitutive_titanmod -#ifndef NEWSTATE - use constitutive_nonlocal -#endif + use constitutive_nonlocal implicit none integer(pInt), parameter :: FILEUNIT = 200_pInt integer(pInt) :: & @@ -176,10 +174,8 @@ subroutine constitutive_init if (any(phase_plasticity == PLASTICITY_J2_ID)) call constitutive_j2_init(FILEUNIT) if (any(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)) call constitutive_phenopowerlaw_init(FILEUNIT) if (any(phase_plasticity == PLASTICITY_DISLOTWIN_ID)) call constitutive_dislotwin_init(FILEUNIT) -#ifndef NEWSTATE if (any(phase_plasticity == PLASTICITY_TITANMOD_ID)) call constitutive_titanmod_init(FILEUNIT) if (any(phase_plasticity == PLASTICITY_NONLOCAL_ID)) call constitutive_nonlocal_init(FILEUNIT) -#endif close(FILEUNIT) write(6,'(/,a)') ' <<<+- constitutive init -+>>>' @@ -431,9 +427,9 @@ subroutine constitutive_init #endif constitutive_sizePostResults(g,i,e) = constitutive_titanmod_sizePostResults(instance) case (PLASTICITY_NONLOCAL_ID) -#ifndef NEWSTATE nonlocalConstitutionPresent = .true. if(myNgrains/=1_pInt) call IO_error(252_pInt, e,i,g) +#ifndef NEWSTATE allocate(constitutive_state0(g,i,e)%p(constitutive_nonlocal_sizeState(instance))) allocate(constitutive_partionedState0(g,i,e)%p(constitutive_nonlocal_sizeState(instance))) allocate(constitutive_subState0(g,i,e)%p(constitutive_nonlocal_sizeState(instance))) @@ -458,20 +454,24 @@ subroutine constitutive_init constitutive_aTolState(g,i,e)%p = constitutive_nonlocal_aTolState(instance) constitutive_sizeState(g,i,e) = constitutive_nonlocal_sizeState(instance) constitutive_sizeDotState(g,i,e) = constitutive_nonlocal_sizeDotState(instance) - constitutive_sizePostResults(g,i,e) = constitutive_nonlocal_sizePostResults(instance) #endif -! constitutive_sizePostResults(g,i,e) = constitutive_nonlocal_sizePostResults(instance) + constitutive_sizePostResults(g,i,e) = constitutive_nonlocal_sizePostResults(instance) end select enddo GrainLoop enddo IPloop enddo ElemLoop -#ifndef NEWSTATE - if (nonlocalConstitutionPresent) & + + if (nonlocalConstitutionPresent) & +#ifdef NEWSTATE + call constitutive_nonlocal_stateInit(mappingConstitutive) +#else call constitutive_nonlocal_stateInit(constitutive_state0(1,1:iMax,1:eMax)) +#endif + +#ifndef NEWSTATE do e = 1_pInt,mesh_NcpElems ! loop over elements myNgrains = homogenization_Ngrains(mesh_element(3,e)) forall(i = 1_pInt:FE_Nips(FE_geomtype(mesh_element(2,e))), g = 1_pInt:myNgrains) - constitutive_partionedState0(g,i,e)%p = constitutive_state0(g,i,e)%p constitutive_state(g,i,e)%p = constitutive_state0(g,i,e)%p ! need to be defined for first call of constitutive_microstructure in crystallite_init endforall @@ -598,14 +598,11 @@ subroutine constitutive_microstructure(temperature, Fe, Fp, ipc, ip, el) PLASTICITY_NONLOCAL_ID use constitutive_titanmod, only: & constitutive_titanmod_microstructure -#ifndef NEWSTATE use constitutive_nonlocal, only: & constitutive_nonlocal_microstructure -#endif use constitutive_dislotwin, only: & constitutive_dislotwin_microstructure - implicit none integer(pInt), intent(in) :: & ipc, & !< grain number @@ -638,8 +635,10 @@ subroutine constitutive_microstructure(temperature, Fe, Fp, ipc, ip, el) call constitutive_titanmod_microstructure(temperature,constitutive_state(ipc,ip,el), & ipc,ip,el) #endif -#ifndef NEWSTATE case (PLASTICITY_NONLOCAL_ID) +#ifdef NEWSTATE + call constitutive_nonlocal_microstructure(mappingConstitutive,Fe,Fp,ipc,ip,el) +#else call constitutive_nonlocal_microstructure(constitutive_state,Fe,Fp,ipc,ip,el) #endif end select @@ -673,10 +672,9 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, temperature, ip constitutive_dislotwin_LpAndItsTangent use constitutive_titanmod, only: & constitutive_titanmod_LpAndItsTangent -#ifndef NEWSTATE use constitutive_nonlocal, only: & constitutive_nonlocal_LpAndItsTangent -#endif + implicit none integer(pInt), intent(in) :: & ipc, & !< grain number @@ -732,8 +730,11 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, temperature, ip call constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, & temperature,constitutive_state(ipc,ip,el),ipc,ip,el) #endif -#ifndef NEWSTATE case (PLASTICITY_NONLOCAL_ID) +#ifdef NEWSTATE + call constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, & + temperature, mappingConstitutive, ipc,ip,el) +#else call constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, & temperature, constitutive_state(ipc,ip,el), ipc,ip,el) #endif @@ -846,10 +847,9 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, Temperature, constitutive_dislotwin_dotState use constitutive_titanmod, only: & constitutive_titanmod_dotState -#ifndef NEWSTATE use constitutive_nonlocal, only: & constitutive_nonlocal_dotState -#endif + implicit none integer(pInt), intent(in) :: & ipc, & !< grain number @@ -914,8 +914,14 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, Temperature, constitutive_dotState(ipc,ip,el)%p = constitutive_titanmod_dotState(Tstar_v,Temperature,& constitutive_state(ipc,ip,el), ipc,ip,el) #endif -#ifndef NEWSTATE case (PLASTICITY_NONLOCAL_ID) +#ifdef NEWSTATE + plasticState(mappingConstitutive(2,ipc,ip,el))%dotState(:,mappingConstitutive(1,ipc,ip,el)) = & + constitutive_nonlocal_dotState(Tstar_v, FeArray, FpArray, & + Temperature, mappingConstitutive, subdt, & + subfracArray, ipc, ip, el) + +#else constitutive_dotState(ipc,ip,el)%p = constitutive_nonlocal_dotState(Tstar_v, FeArray, FpArray, & Temperature, constitutive_state, constitutive_state0, subdt, & subfracArray, ipc, ip, el) @@ -931,10 +937,8 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, Temperature, if (tock < tick) debug_cumDotStateTicks = debug_cumDotStateTicks + maxticks !$OMP END CRITICAL (debugTimingDotState) endif - end subroutine constitutive_collectDotState -#ifndef NEWSTATE !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the incremental change of !> microstructure based on the current stress and state @@ -951,7 +955,10 @@ subroutine constitutive_collectDeltaState(Tstar_v, ipc, ip, el) use material, only: & phase_plasticity, & material_phase, & - PLASTICITY_NONLOCAL_ID +#ifdef NEWSTATE + plasticState, & +#endif + PLASTICITY_NONLOCAL_ID use constitutive_nonlocal, only: & constitutive_nonlocal_deltaState @@ -973,11 +980,17 @@ subroutine constitutive_collectDeltaState(Tstar_v, ipc, ip, el) select case (phase_plasticity(material_phase(ipc,ip,el))) case (PLASTICITY_NONLOCAL_ID) +#ifdef NEWSTATE + call constitutive_nonlocal_deltaState(mappingConstitutive, Tstar_v,ipc,ip,el) +#else call constitutive_nonlocal_deltaState(constitutive_deltaState(ipc,ip,el),& constitutive_state(ipc,ip,el), Tstar_v,ipc,ip,el) - +#endif case default -#ifndef NEWSTATE +#ifdef NEWSTATE + plasticState(mappingConstitutive(2,ipc,ip,el))%deltaState(:,mappingConstitutive(2,ipc,ip,el)) = & + 0.0_pReal +#else constitutive_deltaState(ipc,ip,el)%p = 0.0_pReal !ToDo: needed or will it remain zero anyway? #endif end select @@ -993,8 +1006,6 @@ subroutine constitutive_collectDeltaState(Tstar_v, ipc, ip, el) endif end subroutine constitutive_collectDeltaState -#endif - !-------------------------------------------------------------------------------------------------- !> @brief returns array of constitutive results !-------------------------------------------------------------------------------------------------- @@ -1026,10 +1037,8 @@ function constitutive_postResults(Tstar_v, FeArray, temperature, ipc, ip, el) constitutive_dislotwin_postResults use constitutive_titanmod, only: & constitutive_titanmod_postResults -#ifndef NEWSTATE use constitutive_nonlocal, only: & constitutive_nonlocal_postResults -#endif implicit none integer(pInt), intent(in) :: & ipc, & !< grain number @@ -1084,7 +1093,12 @@ function constitutive_postResults(Tstar_v, FeArray, temperature, ipc, ip, el) #else constitutive_postResults = constitutive_dislotwin_postResults(Tstar_v,Temperature,& constitutive_state(ipc,ip,el),ipc,ip,el) +#endif case (PLASTICITY_NONLOCAL_ID) +#ifdef NEWSTATE + constitutive_postResults = constitutive_nonlocal_postResults(Tstar_v, FeArray, & + mappingConstitutive, ipc, ip, el) +#else constitutive_postResults = constitutive_nonlocal_postResults(Tstar_v, FeArray, & constitutive_state, constitutive_dotstate(ipc,ip,el), ipc, ip, el) #endif diff --git a/code/constitutive_nonlocal.f90 b/code/constitutive_nonlocal.f90 index 7d47fa834..730ea9a54 100644 --- a/code/constitutive_nonlocal.f90 +++ b/code/constitutive_nonlocal.f90 @@ -293,9 +293,20 @@ use material, only: homogenization_maxNgrains, & phase_Noutput, & PLASTICITY_NONLOCAL_label, & PLASTICITY_NONLOCAL_ID, & +#ifdef NEWSTATE + plasticState, & + material_phase, & + material_Nphase, & +#endif MATERIAL_partPhase use lattice +#ifdef NEWSTATE + use numerics,only: & + numerics_integrator +#endif + +implicit none integer(pInt), intent(in) :: fileUnit !*** local variables @@ -322,7 +333,11 @@ integer(pInt) :: phase, & mySize = 0_pInt ! to suppress warnings, safe as init is called only once character(len=65536) :: & tag = '', & - line = '' + line = '' +#ifdef NEWSTATE + integer(pInt) :: sizeState, sizeDotState,sizeDependentState + integer(pInt) :: NofMyPhase +#endif write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_NONLOCAL_label//' init -+>>>' write(6,'(a)') ' $Id$' @@ -838,7 +853,6 @@ allocate(iGamma(maxTotalNslip,maxNinstances), source=0_pInt) allocate(iRhoF(maxTotalNslip,maxNinstances), source=0_pInt) allocate(iTauF(maxTotalNslip,maxNinstances), source=0_pInt) allocate(iTauB(maxTotalNslip,maxNinstances), source=0_pInt) - allocate(burgers(maxTotalNslip,maxNinstances), source=0.0_pReal) allocate(lambda0(maxTotalNslip,maxNinstances), source=0.0_pReal) allocate(minDipoleHeight(maxTotalNslip,2,maxNinstances), source=-1.0_pReal) @@ -870,6 +884,9 @@ allocate(nonSchmidProjection(3,3,4,maxTotalNslip,maxNinstances), initializeInstances: do phase = 1_pInt, size(phase_plasticity) if (phase_plasticity(phase) == PLASTICITY_NONLOCAL_ID) then +#ifdef NEWSTATE + NofMyPhase=count(material_phase==phase) +#endif instance = phase_plasticityInstance(phase) !*** Inverse lookup of my slip system family and the slip system in lattice @@ -890,6 +907,37 @@ allocate(nonSchmidProjection(3,3,4,maxTotalNslip,maxNinstances), constitutive_nonlocal_sizeState(instance) = constitutive_nonlocal_sizeDotState(instance) & + constitutive_nonlocal_sizeDependentState(instance) & + int(size(OTHERSTATES),pInt) * ns + +#ifdef NEWSTATE +! Determine size of state array + plasticState(phase)%nonlocal = .true. + sizeDotState = int(size(BASICSTATES),pInt) * ns + sizeDependentState = int(size(DEPENDENTSTATES),pInt) * ns + sizeState = sizeDotState & + + constitutive_nonlocal_sizeDependentState(instance) & + + int(size(OTHERSTATES),pInt) * ns + + plasticState(phase)%sizeState = sizeState + plasticState(phase)%sizeDotState = sizeDotState + allocate(plasticState(phase)%aTolState (sizeState), source=0.0_pReal) + allocate(plasticState(phase)%state0 (sizeState,NofMyPhase), source=0.0_pReal) + allocate(plasticState(phase)%partionedState0 (sizeState,NofMyPhase), source=0.0_pReal) + allocate(plasticState(phase)%subState0 (sizeState,NofMyPhase), source=0.0_pReal) + 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 (sizeDotState,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) + allocate(plasticState(phase)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal) + endif + if (any(numerics_integrator == 4_pInt)) & + allocate(plasticState(phase)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal) + if (any(numerics_integrator == 5_pInt)) & + allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal) +#endif !*** determine indices to state array @@ -1111,16 +1159,18 @@ allocate(nonSchmidProjection(3,3,4,maxTotalNslip,maxNinstances), nonSchmidProjection(1:3,1:3,t,s,instance) = nonSchmidProjection(1:3,1:3,t,s,instance) & + lattice_Sslip(1:3,1:3,1,slipSystemLattice(s,instance),phase) enddo - endif + endif +#ifdef NEWSTATE + call constitutive_nonlocal_aTolState(phase,instance) +#endif enddo initializeInstances end subroutine constitutive_nonlocal_init - - +#ifdef NEWSTATE !-------------------------------------------------------------------------------------------------- !> @brief sets the initial microstructural state for a given instance of this plasticity !-------------------------------------------------------------------------------------------------- -subroutine constitutive_nonlocal_stateInit(state) +subroutine constitutive_nonlocal_stateInit(mapping) use IO, only: IO_error use lattice, only: lattice_maxNslipFamily use math, only: math_sampleGaussVar @@ -1132,16 +1182,18 @@ use mesh, only: mesh_ipVolume, & FE_geomtype use material, only: material_phase, & phase_plasticityInstance, & + plasticState, & + material_Nphase, & + homogenization_maxNgrains, & phase_plasticity ,& PLASTICITY_NONLOCAL_ID +use numerics,only: & + numerics_integrator implicit none -!*** input/output variables -type(p_vec), dimension(1,mesh_maxNips,mesh_NcpElems), intent(inout) :: & - state ! microstructural state - -!*** local variables +integer(pInt),dimension(2,homogenization_maxngrains,mesh_maxNips,mesh_ncpelems), intent(in) :: & + mapping integer(pInt) el, & ip, & e, & @@ -1162,9 +1214,186 @@ real(pReal) meanDensity, & densityBinning, & minimumIpVolume - maxNinstances = int(count(phase_plasticity == PLASTICITY_NONLOCAL_ID),pInt) +do instance = 1_pInt,maxNinstances + ns = totalNslip(instance) + + ! randomly distribute dislocation segments on random slip system and of random type in the volume + if (rhoSglRandom(instance) > 0.0_pReal) then + + ! get the total volume of the instance + + minimumIpVolume = 1e99_pReal !use huge here? + totalVolume = 0.0_pReal + do e = 1_pInt,mesh_NcpElems + do i = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,e))) + if (PLASTICITY_NONLOCAL_ID == phase_plasticity(material_phase(1,i,e)) & + .and. instance == phase_plasticityInstance(material_phase(1,i,e))) then + totalVolume = totalVolume + mesh_ipVolume(i,e) + minimumIpVolume = min(minimumIpVolume, mesh_ipVolume(i,e)) + endif + enddo + enddo + densityBinning = rhoSglRandomBinning(instance) / minimumIpVolume ** (2.0_pReal / 3.0_pReal) + + ! subsequently fill random ips with dislocation segments until we reach the desired overall density + + meanDensity = 0.0_pReal + do while(meanDensity < rhoSglRandom(instance)) + call random_number(rnd) + el = nint(rnd(1)*real(mesh_NcpElems,pReal)+0.5_pReal,pInt) + ip = nint(rnd(2)*real(FE_Nips(FE_geomtype(mesh_element(2,el))),pReal)+0.5_pReal,pInt) + if (PLASTICITY_NONLOCAL_ID == phase_plasticity(material_phase(1,ip,el)) & + .and. instance == phase_plasticityInstance(material_phase(1,ip,el))) then + s = nint(rnd(3)*real(ns,pReal)+0.5_pReal,pInt) + t = nint(rnd(4)*4.0_pReal+0.5_pReal,pInt) + meanDensity = meanDensity + densityBinning * mesh_ipVolume(ip,el) / totalVolume + plasticState(mapping(2,1,ip,el))%state(iRhoU(s,t,instance),mapping(1,1,ip,el)) = & + plasticState(mapping(2,1,ip,el))%state(iRhoU(s,t,instance),mapping(1,1,ip,el)) & + + densityBinning + + endif + enddo + ! homogeneous distribution of density with some noise + else + do e = 1_pInt,mesh_NcpElems + do i = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,e))) + if (PLASTICITY_NONLOCAL_ID == phase_plasticity(material_phase(1,i,e)) & + .and. instance == phase_plasticityInstance(material_phase(1,i,e))) then + do f = 1_pInt,lattice_maxNslipFamily + from = 1_pInt + sum(Nslip(1:f-1_pInt,instance)) + upto = sum(Nslip(1:f,instance)) + do s = from,upto + do j = 1_pInt,2_pInt + noise(j) = math_sampleGaussVar(0.0_pReal, rhoSglScatter(instance)) + enddo + + plasticState(mapping(2,1,i,e))%state(iRhoU(s,1,instance),mapping(1,1,i,e)) = & + rhoSglEdgePos0(f,instance) + noise(1) + plasticState(mapping(2,1,i,e))%state(iRhoU(s,2,instance),mapping(1,1,i,e)) = & + rhoSglEdgeNeg0(f,instance) + noise(1) + plasticState(mapping(2,1,i,e))%state(iRhoU(s,3,instance),mapping(1,1,i,e)) = & + rhoSglScrewPos0(f,instance) + noise(2) + plasticState(mapping(2,1,i,e))%state(iRhoU(s,4,instance),mapping(1,1,i,e)) = & + rhoSglScrewNeg0(f,instance) + noise(2) + + enddo + + plasticState(mapping(2,1,i,e))%state(iRhoD(from:upto,1,instance),mapping(1,1,i,e)) = & + rhoDipEdge0(f,instance) + plasticState(mapping(2,1,i,e))%state(iRhoD(from:upto,2,instance),mapping(1,1,i,e)) = & + rhoDipScrew0(f,instance) + + enddo + endif + enddo + enddo +do e = 1_pInt,mesh_NcpElems + do i = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,e))) + if (PLASTICITY_NONLOCAL_ID == phase_plasticity(material_phase(1,i,e)) & + .and. instance == phase_plasticityInstance(material_phase(1,i,e))) then + plasticState(mapping(2,1,i,e))%state0(:,mapping(1,1,i,e)) = & + plasticState(mapping(2,1,i,e))%state(:,mapping(1,1,i,e)) + plasticState(mapping(2,1,i,e))%partionedState0(:,mapping(1,1,i,e)) = & + plasticState(mapping(2,1,i,e))%state0(:,mapping(1,1,i,e)) + plasticState(mapping(2,1,i,e))%State(:,mapping(1,1,i,e)) = & + plasticState(mapping(2,1,i,e))%state0(:,mapping(1,1,i,e)) + endif + enddo +enddo + + endif +enddo + +end subroutine constitutive_nonlocal_stateInit + +!-------------------------------------------------------------------------------------------------- +!> @brief sets the relevant state values for a given instance of this plasticity +!-------------------------------------------------------------------------------------------------- +subroutine constitutive_nonlocal_aTolState(phase,instance) +use material, only: & + plasticState +implicit none +!*** input variables +integer(pInt), intent(in) :: & + instance, & !< number specifying the instance of the plasticity + phase + real(pReal), dimension(size(plasticState(phase)%aTolState(:))) :: tempTol +!*** local variables + integer(pInt) :: ns, t, c + + + tempTol = 0.0_pReal + ns = totalNslip(instance) +! constitutive_nonlocal_aTolState = 0.0_pReal +forall (t = 1_pInt:4_pInt) + tempTol(iRhoU(1:ns,t,instance)) = aTolRho(instance) + tempTol(iRhoB(1:ns,t,instance)) = aTolRho(instance) +endforall +forall (c = 1_pInt:2_pInt) & + tempTol(iRhoD(1:ns,c,instance)) = aTolRho(instance) + tempTol(iGamma(1:ns,instance)) = aTolShear(instance) + plasticState(phase)%aTolState = tempTol + +end subroutine constitutive_nonlocal_aTolState + +#else +!-------------------------------------------------------------------------------------------------- +!> @brief sets the initial microstructural state for a given instance of this plasticity +!-------------------------------------------------------------------------------------------------- +subroutine constitutive_nonlocal_stateInit(state) +use IO, only: IO_error +use lattice, only: lattice_maxNslipFamily +use math, only: math_sampleGaussVar +use mesh, only: mesh_ipVolume, & + mesh_NcpElems, & + mesh_maxNips, & + mesh_element, & + FE_Nips, & + FE_geomtype +use material, only: material_phase, & + phase_plasticityInstance, & +#ifdef NEWSTATE + plasticState, & + material_Nphase, & +#endif + phase_plasticity ,& + PLASTICITY_NONLOCAL_ID +#ifdef NEWSTATE + use numerics,only: & + numerics_integrator +#endif + +implicit none + +!*** input/output variables +type(p_vec), dimension(1,mesh_maxNips,mesh_NcpElems), intent(inout) :: & + state ! microstructural state + +!*** local variables +integer(pInt) el, & + ip, & + e, & + i, & + ns, & ! short notation for total number of active slip systems + f, & ! index of lattice family + from, & + upto, & + s, & ! index of slip system + t, & + j, & + instance, & + phase, & + maxNinstances +real(pReal), dimension(2) :: noise +real(pReal), dimension(4) :: rnd +real(pReal) meanDensity, & + totalVolume, & + densityBinning, & + minimumIpVolume + +maxNinstances = int(count(phase_plasticity == PLASTICITY_NONLOCAL_ID),pInt) ! ititalize all states to zero @@ -1212,7 +1441,6 @@ do instance = 1_pInt,maxNinstances state(1,ip,el)%p(iRhoU(s,t,instance)) = state(1,ip,el)%p(iRhoU(s,t,instance)) + densityBinning endif enddo - ! homogeneous distribution of density with some noise else do e = 1_pInt,mesh_NcpElems @@ -1243,7 +1471,6 @@ enddo end subroutine constitutive_nonlocal_stateInit - !-------------------------------------------------------------------------------------------------- !> @brief sets the relevant state values for a given instance of this plasticity !-------------------------------------------------------------------------------------------------- @@ -1272,7 +1499,317 @@ constitutive_nonlocal_aTolState(iGamma(1:ns,instance)) = aTolShear(instance) end function constitutive_nonlocal_aTolState +#endif +#ifdef NEWSTATE +!-------------------------------------------------------------------------------------------------- +!> @brief calculates quantities characterizing the microstructure +!-------------------------------------------------------------------------------------------------- +subroutine constitutive_nonlocal_microstructure(mapping, Fe, Fp, ipc, ip, el) + +use IO, only: & + IO_error +use math, only: & + pi, & + math_mul33x3, & + math_mul3x3, & + math_norm3, & + math_invert33, & + math_transpose33 +use debug, only: & + debug_level, & + debug_constitutive, & + debug_levelBasic, & + debug_levelExtensive, & + debug_levelSelective, & + debug_g, & + debug_i, & + debug_e +use mesh, only: & + mesh_NcpElems, & + mesh_maxNips, & + mesh_element, & + mesh_ipNeighborhood, & + mesh_ipCoordinates, & + mesh_ipVolume, & + mesh_ipAreaNormal, & + mesh_ipArea, & + FE_NipNeighbors, & + mesh_maxNipNeighbors, & + FE_geomtype, & + FE_celltype +use material, only: & + homogenization_maxNgrains, & + material_phase, & + plasticState ,& + phase_localPlasticity, & + phase_plasticityInstance +use lattice, only: & + lattice_sd, & + lattice_st, & + lattice_mu, & + lattice_nu, & + lattice_structure, & + LATTICE_bcc_ID, & + LATTICE_fcc_ID + +implicit none + +!*** input variables +integer(pInt), intent(in) :: ipc, & ! current grain ID + ip, & ! current integration point + el ! current element +real(pReal), dimension(3,3), intent(in) :: & + Fe, & ! elastic deformation gradient + Fp ! elastic deformation gradient + + +!*** output variables + +!*** local variables +integer(pInt) neighbor_el, & ! element number of neighboring material point + neighbor_ip, & ! integration point of neighboring material point + instance, & ! my instance of this plasticity + neighbor_instance, & ! instance of this plasticity of neighboring material point + phase, & + neighbor_phase, & + ns, & ! total number of active slip systems at my material point + neighbor_ns, & ! total number of active slip systems at neighboring material point + c, & ! index of dilsocation character (edge, screw) + s, & ! slip system index + t, & ! index of dilsocation type (e+, e-, s+, s-, used e+, used e-, used s+, used s-) + dir, & + n, & + nRealNeighbors ! number of really existing neighbors +integer(pInt), dimension(2) :: neighbors +real(pReal) detFe, & + detFp, & + FVsize, & + temp, & + correction, & + myRhoForest +real(pReal), dimension(2) :: rhoExcessGradient, & + rhoExcessGradient_over_rho, & + rhoTotal +real(pReal), dimension(3) :: rhoExcessDifferences, & + normal_latticeConf +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & + rhoForest, & ! forest dislocation density + tauBack, & ! back stress from pileup on same slip system + tauThreshold ! threshold shear stress +real(pReal), dimension(3,3) :: invFe, & ! inverse of elastic deformation gradient + invFp, & ! inverse of plastic deformation gradient + connections, & + invConnections +real(pReal), dimension(3,mesh_maxNipNeighbors) :: & + connection_latticeConf +real(pReal), dimension(2,totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & + rhoExcess +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),2) :: & + rhoDip ! dipole dislocation density (edge, screw) +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),8) :: & + rhoSgl ! single dislocation density (edge+, edge-, screw+, screw-, used edge+, used edge-, used screw+, used screw-) +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))), & + totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & + myInteractionMatrix ! corrected slip interaction matrix +real(pReal), dimension(2,maxval(totalNslip),mesh_maxNipNeighbors) :: & + neighbor_rhoExcess, & ! excess density at neighboring material point + neighbor_rhoTotal ! total density at neighboring material point +real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),2) :: & + m ! direction of dislocation motion +logical inversionError +integer(pInt),dimension(2,homogenization_maxngrains,mesh_maxNips,mesh_ncpelems), intent(in) :: & + mapping + +phase = material_phase(ipc,ip,el) +instance = phase_plasticityInstance(phase) +ns = totalNslip(instance) + +!*** get basic states + +forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) + + rhoSgl(s,t) = max(plasticState(mapping(2,1,ip,el))%state(iRhoU(s,t,instance),mapping(1,1,ip,el)), & + 0.0_pReal) ! ensure positive single mobile densities + rhoSgl(s,t+4_pInt) = plasticState(mapping(2,1,ip,el))%state(iRhoB(s,t,instance),mapping(1,1,ip,el)) +endforall +forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) & + rhoDip(s,c) = max(plasticState(mapping(2,1,ip,el))%state(iRhoD(s,c,instance),mapping(1,1,ip,el)), 0.0_pReal) ! ensure positive dipole densities +where (abs(rhoSgl) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(instance) & + .or. abs(rhoSgl) < significantRho(instance)) & + rhoSgl = 0.0_pReal +where (abs(rhoDip) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(instance) & + .or. abs(rhoDip) < significantRho(instance)) & + rhoDip = 0.0_pReal + + +!*** calculate the forest dislocation density +!*** (= projection of screw and edge dislocations) + +forall (s = 1_pInt:ns) & + rhoForest(s) = dot_product((sum(abs(rhoSgl(1:ns,[1,2,5,6])),2) + rhoDip(1:ns,1)), & + forestProjectionEdge(s,1:ns,instance)) & + + dot_product((sum(abs(rhoSgl(1:ns,[3,4,7,8])),2) + rhoDip(1:ns,2)), & + forestProjectionScrew(s,1:ns,instance)) + +!*** calculate the threshold shear stress for dislocation slip +!*** coefficients are corrected for the line tension effect +!*** (see Kubin,Devincre,Hoc; 2008; Modeling dislocation storage rates and mean free paths in face-centered cubic crystals) + +myInteractionMatrix = 0.0_pReal +myInteractionMatrix(1:ns,1:ns) = interactionMatrixSlipSlip(1:ns,1:ns,instance) +if (lattice_structure(phase) == LATTICE_bcc_ID .or. lattice_structure(phase) == LATTICE_fcc_ID) then ! only fcc and bcc + do s = 1_pInt,ns + myRhoForest = max(rhoForest(s),significantRho(instance)) + correction = ( 1.0_pReal - linetensionEffect(instance) & + + linetensionEffect(instance) & + * log(0.35_pReal * burgers(s,instance) * sqrt(myRhoForest)) & + / log(0.35_pReal * burgers(s,instance) * 1e6_pReal)) ** 2.0_pReal + myInteractionMatrix(s,1:ns) = correction * myInteractionMatrix(s,1:ns) + enddo +endif +forall (s = 1_pInt:ns) & + tauThreshold(s) = lattice_mu(phase) * burgers(s,instance) & + * sqrt(dot_product((sum(abs(rhoSgl),2) + sum(abs(rhoDip),2)), myInteractionMatrix(s,1:ns))) + +!*** calculate the dislocation stress of the neighboring excess dislocation densities +!*** zero for material points of local plasticity + +tauBack = 0.0_pReal + +if (.not. phase_localPlasticity(phase) .and. shortRangeStressCorrection(instance)) then + call math_invert33(Fe, invFe, detFe, inversionError) + call math_invert33(Fp, invFp, detFp, inversionError) + rhoExcess(1,1:ns) = rhoSgl(1:ns,1) - rhoSgl(1:ns,2) + rhoExcess(2,1:ns) = rhoSgl(1:ns,3) - rhoSgl(1:ns,4) + FVsize = mesh_ipVolume(ip,el) ** (1.0_pReal/3.0_pReal) + + !* loop through my neighborhood and get the connection vectors (in lattice frame) and the excess densities + + nRealNeighbors = 0_pInt + neighbor_rhoTotal = 0.0_pReal + do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el)))) + neighbor_el = mesh_ipNeighborhood(1,n,ip,el) + neighbor_ip = mesh_ipNeighborhood(2,n,ip,el) + if (neighbor_el > 0 .and. neighbor_ip > 0) then + neighbor_phase = material_phase(ipc,neighbor_ip,neighbor_el) + neighbor_instance = phase_plasticityInstance(neighbor_phase) + neighbor_ns = totalNslip(neighbor_instance) + if (.not. phase_localPlasticity(neighbor_phase) & + .and. neighbor_instance == instance) then ! same instance should be same structure + if (neighbor_ns == ns) then + nRealNeighbors = nRealNeighbors + 1_pInt + forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) + neighbor_rhoExcess(c,s,n) = & + max(plasticState(mapping(2,ipc,neighbor_ip,neighbor_el))% & + state(iRhoU(s,2*c-1,neighbor_instance), & + mapping(1,ipc,neighbor_ip,neighbor_el)),0.0_pReal) &! positive mobiles + - max(plasticState(mapping(2,ipc,neighbor_ip,neighbor_el))% & + state(iRhoU(s,2*c,neighbor_instance), & + mapping(1,ipc,neighbor_ip,neighbor_el)),0.0_pReal) ! negative mobiles + endforall + connection_latticeConf(1:3,n) = & + math_mul33x3(invFe, mesh_ipCoordinates(1:3,neighbor_ip,neighbor_el) & + - mesh_ipCoordinates(1:3,ip,el)) + normal_latticeConf = math_mul33x3(math_transpose33(invFp), mesh_ipAreaNormal(1:3,n,ip,el)) + if (math_mul3x3(normal_latticeConf,connection_latticeConf(1:3,n)) < 0.0_pReal) then ! neighboring connection points in opposite direction to face normal: must be periodic image + connection_latticeConf(1:3,n) = normal_latticeConf * mesh_ipVolume(ip,el) & + / mesh_ipArea(n,ip,el) ! instead take the surface normal scaled with the diameter of the cell + endif + else + ! different number of active slip systems + call IO_error(-1_pInt,ext_msg='different number of active slip systems in neighboring IPs of same crystal structure') + endif + else + ! local neighbor or different lattice structure or different constitution instance -> use central values instead + connection_latticeConf(1:3,n) = 0.0_pReal + neighbor_rhoExcess(1:2,1:ns,n) = rhoExcess + endif + else + ! free surface -> use central values instead + connection_latticeConf(1:3,n) = 0.0_pReal + neighbor_rhoExcess(1:2,1:ns,n) = rhoExcess + endif + enddo + + + !* loop through the slip systems and calculate the dislocation gradient by + !* 1. interpolation of the excess density in the neighorhood + !* 2. interpolation of the dead dislocation density in the central volume + + m(1:3,1:ns,1) = lattice_sd(1:3,slipSystemLattice(1:ns,instance),phase) + m(1:3,1:ns,2) = -lattice_st(1:3,slipSystemLattice(1:ns,instance),phase) + + do s = 1_pInt,ns + + !* gradient from interpolation of neighboring excess density + + do c = 1_pInt,2_pInt + do dir = 1_pInt,3_pInt + neighbors(1) = 2_pInt * dir - 1_pInt + neighbors(2) = 2_pInt * dir + connections(dir,1:3) = connection_latticeConf(1:3,neighbors(1)) & + - connection_latticeConf(1:3,neighbors(2)) + rhoExcessDifferences(dir) = neighbor_rhoExcess(c,s,neighbors(1)) & + - neighbor_rhoExcess(c,s,neighbors(2)) + enddo + call math_invert33(connections,invConnections,temp,inversionError) + if (inversionError) then + call IO_error(-1_pInt,ext_msg='back stress calculation: inversion error') + endif + rhoExcessGradient(c) = math_mul3x3(m(1:3,s,c), & + math_mul33x3(invConnections,rhoExcessDifferences)) + enddo + + !* plus gradient from deads + + do t = 1_pInt,4_pInt + c = (t - 1_pInt) / 2_pInt + 1_pInt + rhoExcessGradient(c) = rhoExcessGradient(c) + rhoSgl(s,t+4_pInt) / FVsize + enddo + + !* normalized with the total density + + rhoExcessGradient_over_rho = 0.0_pReal + forall (c = 1_pInt:2_pInt) & + rhoTotal(c) = (sum(abs(rhoSgl(s,[2*c-1,2*c,2*c+3,2*c+4]))) + rhoDip(s,c) & + + sum(neighbor_rhoTotal(c,s,:))) / real(1_pInt + nRealNeighbors,pReal) + forall (c = 1_pInt:2_pInt, rhoTotal(c) > 0.0_pReal) & + rhoExcessGradient_over_rho(c) = rhoExcessGradient(c) / rhoTotal(c) + + !* gives the local stress correction when multiplied with a factor + + tauBack(s) = - lattice_mu(phase) * burgers(s,instance) / (2.0_pReal * pi) & + * (rhoExcessGradient_over_rho(1) / (1.0_pReal - lattice_nu(phase)) & + + rhoExcessGradient_over_rho(2)) + + enddo +endif + + +!*** set dependent states + +plasticState(mapping(2,ipc,ip,el))%state(iRhoF(1:ns,instance),mapping(1,ipc,ip,el)) = rhoForest +plasticState(mapping(2,ipc,ip,el))%state(iTauF(1:ns,instance),mapping(1,ipc,ip,el)) = tauThreshold +plasticState(mapping(2,ipc,ip,el))%state(iTauB(1:ns,instance),mapping(1,ipc,ip,el)) = tauBack + +#ifndef _OPENMP + if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt & + .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)& + .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then + write(6,*) + write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_microstructure at el ip g',el,ip,ipc + write(6,*) + write(6,'(a,/,12x,12(e10.3,1x))') '<< CONST >> rhoForest', rhoForest + write(6,'(a,/,12x,12(f10.5,1x))') '<< CONST >> tauThreshold / MPa', tauThreshold/1e6 + write(6,'(a,/,12x,12(f10.5,1x))') '<< CONST >> tauBack / MPa', tauBack/1e6 + write(6,*) + endif +#endif + +end subroutine constitutive_nonlocal_microstructure + +#else !-------------------------------------------------------------------------------------------------- !> @brief calculates quantities characterizing the microstructure !-------------------------------------------------------------------------------------------------- @@ -1392,12 +1929,10 @@ real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(ipc, m ! direction of dislocation motion logical inversionError - phase = material_phase(ipc,ip,el) instance = phase_plasticityInstance(phase) ns = totalNslip(instance) - !*** get basic states forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) @@ -1424,7 +1959,6 @@ forall (s = 1_pInt:ns) & forestProjectionScrew(s,1:ns,instance)) - !*** calculate the threshold shear stress for dislocation slip !*** coefficients are corrected for the line tension effect !*** (see Kubin,Devincre,Hoc; 2008; Modeling dislocation storage rates and mean free paths in face-centered cubic crystals) @@ -1446,7 +1980,6 @@ forall (s = 1_pInt:ns) & * sqrt(dot_product((sum(abs(rhoSgl),2) + sum(abs(rhoDip),2)), myInteractionMatrix(s,1:ns))) - !*** calculate the dislocation stress of the neighboring excess dislocation densities !*** zero for material points of local plasticity @@ -1587,6 +2120,7 @@ state(ipc,ip,el)%p(iTauB(1:ns,instance)) = tauBack end subroutine constitutive_nonlocal_microstructure +#endif !-------------------------------------------------------------------------------------------------- !> @brief calculates kinetics @@ -1750,7 +2284,212 @@ endif end subroutine constitutive_nonlocal_kinetics +#ifdef NEWSTATE +!-------------------------------------------------------------------------------------------------- +!> @brief calculates plastic velocity gradient and its tangent +!-------------------------------------------------------------------------------------------------- +subroutine constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar99, Tstar_v, Temperature, mapping, ipc, ip, el) +use math, only: math_Plain3333to99, & + math_mul6x6, & + math_mul33xx33, & + math_Mandel6to33 +use debug, only: debug_level, & + debug_constitutive, & + debug_levelBasic, & + debug_levelExtensive, & + debug_levelSelective, & + debug_g, & + debug_i, & + debug_e +use material, only: material_phase, & + phase_plasticityInstance, & + plasticState, & + homogenization_maxngrains +use lattice, only: lattice_Sslip, & + lattice_Sslip_v, & + lattice_NnonSchmid +use mesh, only: mesh_ipVolume, & + mesh_maxNips, & + mesh_ncpelems + +implicit none + +!*** input variables +integer(pInt), intent(in) :: ipc, & !< current grain number + ip, & !< current integration point + el !< current element number +real(pReal), intent(in) :: Temperature !< temperature +real(pReal), dimension(6), intent(in) :: Tstar_v !< 2nd Piola-Kirchhoff stress in Mandel notation + +!*** input/output variables + +!*** output variables +real(pReal), dimension(3,3), intent(out) :: Lp !< plastic velocity gradient +real(pReal), dimension(9,9), intent(out) :: dLp_dTstar99 !< derivative of Lp with respect to Tstar (9x9 matrix) + +!*** local variables +integer(pInt) instance, & !< current instance of this plasticity + phase, & !< phase + ns, & !< short notation for the total number of active slip systems + i, & + j, & + k, & + l, & + t, & !< dislocation type + s, & !< index of my current slip system + sLattice !< index of my current slip system according to lattice order +real(pReal), dimension(3,3,3,3) :: dLp_dTstar3333 !< derivative of Lp with respect to Tstar (3x3x3x3 matrix) +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),8) :: & + rhoSgl !< single dislocation densities (including blocked) +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),4) :: & + v, & !< velocity + tauNS, & !< resolved shear stress including non Schmid and backstress terms + dv_dtau, & !< velocity derivative with respect to the shear stress + dv_dtauNS !< velocity derivative with respect to the shear stress +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & + tau, & !< resolved shear stress including backstress terms + gdotTotal, & !< shear rate + tauBack, & !< back stress from dislocation gradients on same slip system + tauThreshold !< threshold shear stress + +integer(pInt),dimension(2,homogenization_maxngrains,mesh_maxNips,mesh_ncpelems), intent(in) :: & + mapping + +!*** initialize local variables + +Lp = 0.0_pReal +dLp_dTstar3333 = 0.0_pReal + +phase = material_phase(ipc,ip,el) +instance = phase_plasticityInstance(phase) +ns = totalNslip(instance) + + +!*** shortcut to state variables + + +forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) + rhoSgl(s,t) = max(plasticState(mapping(2,1,ip,el))%state(iRhoU(s,t,instance),mapping(1,ipc,ip,el)), & + 0.0_pReal) ! ensure positive single mobile densities + rhoSgl(s,t+4_pInt) = plasticState(mapping(2,1,ip,el))%state(iRhoB(s,t,instance),mapping(1,ipc,ip,el)) +endforall +where (abs(rhoSgl) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(instance) & + .or. abs(rhoSgl) < significantRho(instance)) & + rhoSgl = 0.0_pReal + +tauBack = plasticState(mapping(2,1,ip,el))%state(iTauB(1:ns,instance),mapping(1,ipc,ip,el)) +tauThreshold = plasticState(mapping(2,1,ip,el))%state(iTauF(1:ns,instance),mapping(1,ipc,ip,el)) + + +!*** get resolved shear stress +!*** for screws possible non-schmid contributions are also taken into account + +do s = 1_pInt,ns + sLattice = slipSystemLattice(s,instance) + tau(s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,phase)) + tauNS(s,1) = tau(s) + tauNS(s,2) = tau(s) + if (tau(s) > 0.0_pReal) then + tauNS(s,3) = math_mul33xx33(math_Mandel6to33(Tstar_v), nonSchmidProjection(1:3,1:3,1,s,instance)) + tauNS(s,4) = math_mul33xx33(math_Mandel6to33(Tstar_v), nonSchmidProjection(1:3,1:3,3,s,instance)) + else + tauNS(s,3) = math_mul33xx33(math_Mandel6to33(Tstar_v), nonSchmidProjection(1:3,1:3,2,s,instance)) + tauNS(s,4) = math_mul33xx33(math_Mandel6to33(Tstar_v), nonSchmidProjection(1:3,1:3,4,s,instance)) + endif +enddo +forall (t = 1_pInt:4_pInt) & + tauNS(1:ns,t) = tauNS(1:ns,t) + tauBack ! add backstress +tau = tau + tauBack ! add backstress + + +!*** get dislocation velocity and its tangent and store the velocity in the state array + +! edges +call constitutive_nonlocal_kinetics(v(1:ns,1), dv_dtau(1:ns,1), dv_dtauNS(1:ns,1), & + tau(1:ns), tauNS(1:ns,1), tauThreshold(1:ns), & + 1_pInt, Temperature, ipc, ip, el) +v(1:ns,2) = v(1:ns,1) +dv_dtau(1:ns,2) = dv_dtau(1:ns,1) +dv_dtauNS(1:ns,2) = dv_dtauNS(1:ns,1) + +!screws +if (lattice_NnonSchmid(phase) == 0_pInt) then ! no non-Schmid contributions + forall(t = 3_pInt:4_pInt) + v(1:ns,t) = v(1:ns,1) + dv_dtau(1:ns,t) = dv_dtau(1:ns,1) + dv_dtauNS(1:ns,t) = dv_dtauNS(1:ns,1) + endforall +else ! take non-Schmid contributions into account + do t = 3_pInt,4_pInt + call constitutive_nonlocal_kinetics(v(1:ns,t), dv_dtau(1:ns,t), dv_dtauNS(1:ns,t), & + tau(1:ns), tauNS(1:ns,t), tauThreshold(1:ns), & + 2_pInt , Temperature, ipc, ip, el) + enddo +endif + + +!*** store velocity in state + +forall (t = 1_pInt:4_pInt) & + plasticState(mapping(2,ipc,ip,el))%state(iV(1:ns,t,instance),mapping(1,ipc,ip,el)) = v(1:ns,t) + + +!*** Bauschinger effect + +forall (s = 1_pInt:ns, t = 5_pInt:8_pInt, rhoSgl(s,t) * v(s,t-4_pInt) < 0.0_pReal) & + rhoSgl(s,t-4_pInt) = rhoSgl(s,t-4_pInt) + abs(rhoSgl(s,t)) + + +!*** Calculation of Lp and its tangent + +gdotTotal = sum(rhoSgl(1:ns,1:4) * v, 2) * burgers(1:ns,instance) + +do s = 1_pInt,ns + sLattice = slipSystemLattice(s,instance) + Lp = Lp + gdotTotal(s) * lattice_Sslip(1:3,1:3,1,sLattice,phase) + + ! Schmid contributions to tangent + forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt,k=1_pInt:3_pInt,l=1_pInt:3_pInt) & + dLp_dTstar3333(i,j,k,l) = dLp_dTstar3333(i,j,k,l) & + + lattice_Sslip(i,j,1,sLattice,phase) * lattice_Sslip(k,l,1,sLattice,phase) & + * sum(rhoSgl(s,1:4) * dv_dtau(s,1:4)) * burgers(s,instance) + + ! non Schmid contributions to tangent + if (tau(s) > 0.0_pReal) then + forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt,k=1_pInt:3_pInt,l=1_pInt:3_pInt) & + dLp_dTstar3333(i,j,k,l) = dLp_dTstar3333(i,j,k,l) & + + lattice_Sslip(i,j,1,sLattice,phase) & + * ( nonSchmidProjection(k,l,1,s,instance) * rhoSgl(s,3) * dv_dtauNS(s,3) & + + nonSchmidProjection(k,l,3,s,instance) * rhoSgl(s,4) * dv_dtauNS(s,4) ) & + * burgers(s,instance) + else + forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt,k=1_pInt:3_pInt,l=1_pInt:3_pInt) & + dLp_dTstar3333(i,j,k,l) = dLp_dTstar3333(i,j,k,l) & + + lattice_Sslip(i,j,1,sLattice,phase) & + * ( nonSchmidProjection(k,l,2,s,instance) * rhoSgl(s,3) * dv_dtauNS(s,3) & + + nonSchmidProjection(k,l,4,s,instance) * rhoSgl(s,4) * dv_dtauNS(s,4) ) & + * burgers(s,instance) + endif +enddo +dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333) + + +#ifndef _OPENMP + if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt & + .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)& + .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then + write(6,*) + write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_LpandItsTangent at el ip ipc ',el,ip,ipc + write(6,*) + write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> gdot total / 1e-3',gdotTotal*1e3_pReal + write(6,'(a,/,3(12x,3(f12.7,1x),/))') '<< CONST >> Lp',transpose(Lp) + endif +#endif + +end subroutine constitutive_nonlocal_LpAndItsTangent + +#else !-------------------------------------------------------------------------------------------------- !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- @@ -1948,8 +2687,199 @@ dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333) #endif end subroutine constitutive_nonlocal_LpAndItsTangent +#endif + +#ifdef NEWSTATE +!-------------------------------------------------------------------------------------------------- +!> @brief (instantaneous) incremental change of microstructure +!-------------------------------------------------------------------------------------------------- +subroutine constitutive_nonlocal_deltaState(mapping, Tstar_v, ipc,ip,el) + +use debug, only: debug_level, & + debug_constitutive, & + debug_levelBasic, & + debug_levelExtensive, & + debug_levelSelective, & + debug_g, & + debug_i, & + debug_e +use math, only: pi, & + math_mul6x6 +use lattice, only: lattice_Sslip_v ,& + lattice_mu, & + lattice_nu +use mesh, only: mesh_NcpElems, & + mesh_maxNips, & + mesh_ipVolume +use material, only: homogenization_maxNgrains, & + material_phase, & + phase_plasticityInstance, & + plasticState + +implicit none + +!*** input variables +integer(pInt), intent(in) :: ipc, & ! current grain number + ip, & ! current integration point + el ! current element number +real(pReal), dimension(6), intent(in) :: Tstar_v ! current 2nd Piola-Kirchhoff stress in Mandel notation + +!*** local variables +integer(pInt) phase, & + instance, & ! current instance of this plasticity + ns, & ! short notation for the total number of active slip systems + c, & ! character of dislocation + t, & ! type of dislocation + s, & ! index of my current slip system + sLattice ! index of my current slip system according to lattice order +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),10) :: & + deltaRho, & ! density increment + deltaRhoRemobilization, & ! density increment by remobilization + deltaRhoDipole2SingleStress ! density increment by dipole dissociation (by stress change) +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),8) :: & + rhoSgl ! current single dislocation densities (positive/negative screw and edge without dipoles) +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),4) :: & + v ! dislocation glide velocity +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & + tau, & ! current resolved shear stress + tauBack ! current back stress from pileups on same slip system +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),2) :: & + rhoDip, & ! current dipole dislocation densities (screw and edge dipoles) + dLower, & ! minimum stable dipole distance for edges and screws + dUpper, & ! current maximum stable dipole distance for edges and screws + dUpperOld, & ! old maximum stable dipole distance for edges and screws + deltaDUpper ! change in maximum stable dipole distance for edges and screws +integer(pInt),dimension(2,homogenization_maxngrains,mesh_maxNips,mesh_ncpelems), intent(in) :: & + mapping +#ifndef _OPENMP + if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt & + .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)& + .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then + write(6,*) + write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_deltaState at el ip ipc ',el,ip,ipc + write(6,*) + endif +#endif + +phase = material_phase(ipc,ip,el) +instance = phase_plasticityInstance(phase) +ns = totalNslip(instance) +!*** shortcut to state variables + +forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) + rhoSgl(s,t) = max(plasticState(mapping(2,ipc,ip,el))%state(iRhoU(s,t,instance), & + mapping(1,ipc,ip,el)),0.0_pReal) ! ensure positive single mobile densities + rhoSgl(s,t+4_pInt) = plasticState(mapping(2,ipc,ip,el))%state(iRhoB(s,t,instance), & + mapping(1,ipc,ip,el)) + v(s,t) = plasticState(mapping(2,ipc,ip,el))%state(iRhoU(s,t,instance),mapping(1,ipc,ip,el)) +endforall +forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) + rhoDip(s,c) = max(plasticState(mapping(2,ipc,ip,el))%state(iRhoD(s,c,instance),mapping(1,ipc,ip,el)), 0.0_pReal) ! ensure positive dipole densities + dUpperOld(s,c) = plasticState(mapping(2,ipc,ip,el))%state(iD(s,c,instance),mapping(1,ipc,ip,el)) + +endforall +tauBack = plasticState(mapping(2,ipc,ip,el))%state(iTauB(1:ns,instance),mapping(1,ipc,ip,el)) + +where (abs(rhoSgl) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(instance) & + .or. abs(rhoSgl) < significantRho(instance)) & + rhoSgl = 0.0_pReal +where (abs(rhoDip) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(instance) & + .or. abs(rhoDip) < significantRho(instance)) & + rhoDip = 0.0_pReal + + +!**************************************************************************** +!*** dislocation remobilization (bauschinger effect) + +deltaRhoRemobilization = 0.0_pReal +do t = 1_pInt,4_pInt + do s = 1_pInt,ns + if (rhoSgl(s,t+4_pInt) * v(s,t) < 0.0_pReal) then + deltaRhoRemobilization(s,t) = abs(rhoSgl(s,t+4_pInt)) + rhoSgl(s,t) = rhoSgl(s,t) + abs(rhoSgl(s,t+4_pInt)) + deltaRhoRemobilization(s,t+4_pInt) = - rhoSgl(s,t+4_pInt) + rhoSgl(s,t+4_pInt) = 0.0_pReal + endif + enddo +enddo + + +!**************************************************************************** +!*** calculate dipole formation and dissociation by stress change + +!*** calculate limits for stable dipole height + +do s = 1_pInt,ns + sLattice = slipSystemLattice(s,instance) + tau(s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,phase)) + tauBack(s) + if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal +enddo +dLower = minDipoleHeight(1:ns,1:2,instance) +dUpper(1:ns,1) = lattice_mu(phase) * burgers(1:ns,instance) & + / (8.0_pReal * pi * (1.0_pReal - lattice_nu(phase)) * abs(tau)) +dUpper(1:ns,2) = lattice_mu(phase) * burgers(1:ns,instance) / (4.0_pReal * pi * abs(tau)) + + +!in the test there is an FPE exception. Division by zero? +forall (c = 1_pInt:2_pInt) & + dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) & + + abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), & + dUpper(1:ns,c)) +dUpper = max(dUpper,dLower) +deltaDUpper = dUpper - dUpperOld + +!*** dissociation by stress increase + +deltaRhoDipole2SingleStress = 0.0_pReal +forall (c=1_pInt:2_pInt, s=1_pInt:ns, deltaDUpper(s,c) < 0.0_pReal) & + deltaRhoDipole2SingleStress(s,8_pInt+c) = rhoDip(s,c) * deltaDUpper(s,c) & + / (dUpperOld(s,c) - dLower(s,c)) + +forall (t=1_pInt:4_pInt) & + deltaRhoDipole2SingleStress(1_pInt:ns,t) = -0.5_pReal & + * deltaRhoDipole2SingleStress(1_pInt:ns,(t-1_pInt)/2_pInt+9_pInt) + + + +!*** store new maximum dipole height in state + +forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) & + + plasticState(mapping(2,ipc,ip,el))%state(iD(s,c,instance),mapping(1,ipc,ip,el)) = dUpper(s,c) + +!**************************************************************************** +!*** assign the changes in the dislocation densities to deltaState + +deltaRho = deltaRhoRemobilization & + + deltaRhoDipole2SingleStress + +!deltaState%p = 0.0_pReal +plasticState(mapping(2,ipc,ip,el))%deltaState(:,mapping(1,ipc,ip,el)) = 0.0_pReal +forall (s = 1:ns, t = 1_pInt:4_pInt) + plasticState(mapping(2,ipc,ip,el))%deltaState(iRhoU(s,t,instance),mapping(1,ipc,ip,el)) = & + deltaRho(s,t) + plasticState(mapping(2,ipc,ip,el))%deltaState(iRhoB(s,t,instance),mapping(1,ipc,ip,el)) = & + deltaRho(s,t+4_pInt) +endforall +forall (s = 1:ns, c = 1_pInt:2_pInt) & + plasticState(mapping(2,ipc,ip,el))%deltaState(iRhoD(s,c,instance),mapping(1,ipc,ip,el)) = & + deltaRho(s,c+8_pInt) + +#ifndef _OPENMP + if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt & + .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)& + .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then + write(6,'(a,/,8(12x,12(e12.5,1x),/))') '<< CONST >> dislocation remobilization', deltaRhoRemobilization(1:ns,1:8) + write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> dipole dissociation by stress increase', deltaRhoDipole2SingleStress + write(6,*) + endif +#endif + +end subroutine constitutive_nonlocal_deltaState + +#else !-------------------------------------------------------------------------------------------------- !> @brief (instantaneous) incremental change of microstructure !-------------------------------------------------------------------------------------------------- @@ -2144,8 +3074,630 @@ forall (s = 1:ns, c = 1_pInt:2_pInt) & end subroutine constitutive_nonlocal_deltaState +#endif + +#ifdef NEWSTATE +!-------------------------------------------------------------------------------------------------- +!> @brief calculates the rate of change of microstructure +!-------------------------------------------------------------------------------------------------- +function constitutive_nonlocal_dotState(Tstar_v, Fe, Fp, Temperature, mapping, timestep, subfrac, ipc,ip,el) + +use prec, only: DAMASK_NaN +use numerics, only: numerics_integrationMode, & + numerics_timeSyncing +use IO, only: IO_error +use debug, only: debug_level, & + debug_constitutive, & + debug_levelBasic, & + debug_levelExtensive, & + debug_levelSelective, & + debug_g, & + debug_i, & + debug_e +use math, only: math_norm3, & + math_mul6x6, & + math_mul3x3, & + math_mul33x3, & + math_mul33x33, & + math_inv33, & + math_det33, & + math_transpose33, & + pi +use mesh, only: mesh_NcpElems, & + mesh_maxNips, & + mesh_element, & + mesh_ipNeighborhood, & + mesh_ipVolume, & + mesh_ipArea, & + mesh_ipAreaNormal, & + FE_NipNeighbors, & + FE_geomtype, & + FE_celltype +use material, only: homogenization_maxNgrains, & + material_phase, & + phase_plasticityInstance, & + phase_localPlasticity, & + phase_plasticity ,& + PLASTICITY_NONLOCAL_ID, & + plasticState +use lattice, only: lattice_Sslip_v, & + lattice_sd, & + lattice_st ,& + lattice_mu, & + lattice_nu, & + lattice_structure, & + LATTICE_bcc_ID, & + LATTICE_fcc_ID + +implicit none + +!*** input variables +integer(pInt), intent(in) :: ipc, & !< current grain number + ip, & !< current integration point + el !< current element number +real(pReal), intent(in) :: Temperature, & !< temperature + timestep !< substepped crystallite time increment +real(pReal), dimension(6), intent(in) :: Tstar_v !< current 2nd Piola-Kirchhoff stress in Mandel notation +real(pReal), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + subfrac !< fraction of timestep at the beginning of the substepped crystallite time increment +real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + Fe, & !< elastic deformation gradient + Fp !< plastic deformation gradient +integer(pInt),dimension(2,homogenization_maxngrains,mesh_maxNips,mesh_ncpelems), intent(in) :: & + mapping +real(pReal), allocatable, dimension(:) :: & + constitutive_nonlocal_dotState !< evolution of state variables / microstructure + +!*** local variables +integer(pInt) :: phase, & + instance, & !< current instance of this plasticity + neighbor_instance, & !< instance of my neighbor's plasticity + ns, & !< short notation for the total number of active slip systems + c, & !< character of dislocation + n, & !< index of my current neighbor + neighbor_el, & !< element number of my neighbor + neighbor_ip, & !< integration point of my neighbor + neighbor_n, & !< neighbor index pointing to me when looking from my neighbor + opposite_neighbor, & !< index of my opposite neighbor + opposite_ip, & !< ip of my opposite neighbor + opposite_el, & !< element index of my opposite neighbor + opposite_n, & !< neighbor index pointing to me when looking from my opposite neighbor + t, & !< type of dislocation + topp, & !< type of dislocation with opposite sign to t + s, & !< index of my current slip system + sLattice !< index of my current slip system according to lattice order +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),10) :: & + rhoDot, & !< density evolution + rhoDotMultiplication, & !< density evolution by multiplication + rhoDotFlux, & !< density evolution by flux + rhoDotSingle2DipoleGlide, & !< density evolution by dipole formation (by glide) + rhoDotAthermalAnnihilation, & !< density evolution by athermal annihilation + rhoDotThermalAnnihilation !< density evolution by thermal annihilation +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),8) :: & + rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles) + rhoSglOriginal, & + neighbor_rhoSgl, & !< current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles) + rhoSgl0, & !< single dislocation densities at start of cryst inc (positive/negative screw and edge without dipoles) + my_rhoSgl !< single dislocation densities of central ip (positive/negative screw and edge without dipoles) +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),4) :: & + v, & !< current dislocation glide velocity + v0, & !< dislocation glide velocity at start of cryst inc + my_v, & !< dislocation glide velocity of central ip + neighbor_v, & !< dislocation glide velocity of enighboring ip + gdot !< shear rates +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & + rhoForest, & !< forest dislocation density + tauThreshold, & !< threshold shear stress + tau, & !< current resolved shear stress + tauBack, & !< current back stress from pileups on same slip system + vClimb, & !< climb velocity of edge dipoles + nSources +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),2) :: & + rhoDip, & !< current dipole dislocation densities (screw and edge dipoles) + rhoDipOriginal, & + dLower, & !< minimum stable dipole distance for edges and screws + dUpper !< current maximum stable dipole distance for edges and screws +real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),4) :: & + m !< direction of dislocation motion +real(pReal), dimension(3,3) :: my_F, & !< my total deformation gradient + neighbor_F, & !< total deformation gradient of my neighbor + my_Fe, & !< my elastic deformation gradient + neighbor_Fe, & !< elastic deformation gradient of my neighbor + Favg !< average total deformation gradient of me and my neighbor +real(pReal), dimension(3) :: normal_neighbor2me, & !< interface normal pointing from my neighbor to me in neighbor's lattice configuration + normal_neighbor2me_defConf, & !< interface normal pointing from my neighbor to me in shared deformed configuration + normal_me2neighbor, & !< interface normal pointing from me to my neighbor in my lattice configuration + normal_me2neighbor_defConf !< interface normal pointing from me to my neighbor in shared deformed configuration +real(pReal) area, & !< area of the current interface + transmissivity, & !< overall transmissivity of dislocation flux to neighboring material point + lineLength, & !< dislocation line length leaving the current interface + selfDiffusion, & !< self diffusion + rnd, & + meshlength +logical considerEnteringFlux, & + considerLeavingFlux + +#ifndef _OPENMP + if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt & + .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)& + .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then + write(6,*) + write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_dotState at el ip ipc ',el,ip,ipc + write(6,*) + endif +#endif + +phase = material_phase(ipc,ip,el) +instance = phase_plasticityInstance(phase) +ns = totalNslip(instance) + +tau = 0.0_pReal +gdot = 0.0_pReal + +allocate(constitutive_nonlocal_dotState(plasticState(mapping(2,ipc,ip,el))%sizeDotState),source=0.0_pReal) +!*** shortcut to state variables +forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) + rhoSgl(s,t) = max(plasticState(mapping(2,ipc,ip,el))%state(iRhoU(s,t,instance),mapping(1,ipc,ip,el)), 0.0_pReal) ! ensure positive single mobile densities + rhoSgl(s,t+4_pInt) = plasticState(mapping(2,ipc,ip,el))%state(iRhoB(s,t,instance),mapping(1,ipc,ip,el)) + v(s,t) = plasticState(mapping(2,ipc,ip,el))%state(iV(s,t,instance),mapping(1,ipc,ip,el)) +endforall +forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) + rhoDip(s,c) = max(plasticState(mapping(2,ipc,ip,el))%state(iRhoD(s,c,instance),mapping(1,ipc,ip,el)), 0.0_pReal) ! ensure positive dipole densities + +endforall +rhoForest = plasticState(mapping(2,ipc,ip,el))%state(iRhoF(1:ns,instance),mapping(1,ipc,ip,el)) +tauThreshold = plasticState(mapping(2,ipc,ip,el))%state(iTauF(1:ns,instance),mapping(1,ipc,ip,el)) +tauBack = plasticState(mapping(2,ipc,ip,el))%state(iTauB(1:ns,instance),mapping(1,ipc,ip,el)) + +rhoSglOriginal = rhoSgl +rhoDipOriginal = rhoDip +where (abs(rhoSgl) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(instance) & + .or. abs(rhoSgl) < significantRho(instance)) & + rhoSgl = 0.0_pReal +where (abs(rhoDip) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(instance) & + .or. abs(rhoDip) < significantRho(instance)) & + rhoDip = 0.0_pReal + +if (numerics_timeSyncing) then + forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) + + rhoSgl0(s,t) = max(plasticState(mapping(2,ipc,ip,el))%state(iRhoU(s,t,instance), & + mapping(1,ipc,ip,el)), 0.0_pReal) + rhoSgl0(s,t+4_pInt) = plasticState(mapping(2,ipc,ip,el))%state(iRhoB(s,t,instance),mapping(1,ipc,ip,el)) + v0(s,t) = plasticState(mapping(2,ipc,ip,el))%state(iV(s,t,instance),mapping(1,ipc,ip,el)) + endforall + where (abs(rhoSgl0) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(instance) & + .or. abs(rhoSgl0) < significantRho(instance)) & + rhoSgl0 = 0.0_pReal +endif + + + +!*** sanity check for timestep + +if (timestep <= 0.0_pReal) then ! if illegal timestep... + constitutive_nonlocal_dotState = 0.0_pReal ! ...return without doing anything (-> zero dotState) + return +endif + + + +!**************************************************************************** +!*** Calculate shear rate + +forall (t = 1_pInt:4_pInt) & + gdot(1_pInt:ns,t) = rhoSgl(1_pInt:ns,t) * burgers(1:ns,instance) * v(1:ns,t) + +#ifndef _OPENMP + if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt & + .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)& + .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then + write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> rho / 1/m^2', rhoSgl, rhoDip + write(6,'(a,/,4(12x,12(e12.5,1x),/))') '<< CONST >> gdot / 1/s',gdot + endif +#endif + + + +!**************************************************************************** +!*** calculate limits for stable dipole height + +do s = 1_pInt,ns ! loop over slip systems + sLattice = slipSystemLattice(s,instance) + tau(s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,phase)) + tauBack(s) + if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal +enddo + +dLower = minDipoleHeight(1:ns,1:2,instance) +dUpper(1:ns,1) = lattice_mu(phase) * burgers(1:ns,instance) & + / (8.0_pReal * pi * (1.0_pReal - lattice_nu(phase)) * abs(tau)) +dUpper(1:ns,2) = lattice_mu(phase) * burgers(1:ns,instance) & + / (4.0_pReal * pi * abs(tau)) +forall (c = 1_pInt:2_pInt) & + dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) & + + abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), & + dUpper(1:ns,c)) +dUpper = max(dUpper,dLower) + + + +!**************************************************************************** +!*** calculate dislocation multiplication + +rhoDotMultiplication = 0.0_pReal +if (lattice_structure(phase) == LATTICE_bcc_ID) then ! BCC + forall (s = 1:ns, sum(abs(v(s,1:4))) > 0.0_pReal) + rhoDotMultiplication(s,1:2) = sum(abs(gdot(s,3:4))) / burgers(s,instance) & ! assuming double-cross-slip of screws to be decisive for multiplication + * sqrt(rhoForest(s)) / lambda0(s,instance) ! & ! mean free path + ! * 2.0_pReal * sum(abs(v(s,3:4))) / sum(abs(v(s,1:4))) ! ratio of screw to overall velocity determines edge generation + rhoDotMultiplication(s,3:4) = sum(abs(gdot(s,3:4))) / burgers(s,instance) & ! assuming double-cross-slip of screws to be decisive for multiplication + * sqrt(rhoForest(s)) / lambda0(s,instance) ! & ! mean free path + ! * 2.0_pReal * sum(abs(v(s,1:2))) / sum(abs(v(s,1:4))) ! ratio of edge to overall velocity determines screw generation + endforall + +else ! ALL OTHER STRUCTURES + if (probabilisticMultiplication(instance)) then + meshlength = mesh_ipVolume(ip,el)**0.333_pReal + where(sum(rhoSgl(1:ns,1:4),2) > 0.0_pReal) + nSources = (sum(rhoSgl(1:ns,1:2),2) * fEdgeMultiplication(instance) + sum(rhoSgl(1:ns,3:4),2)) & + / sum(rhoSgl(1:ns,1:4),2) * meshlength / lambda0(1:ns,instance) * sqrt(rhoForest(1:ns)) + elsewhere + nSources = meshlength / lambda0(1:ns,instance) * sqrt(rhoForest(1:ns)) + endwhere + do s = 1_pInt,ns + if (nSources(s) < 1.0_pReal) then + if (sourceProbability(s,ipc,ip,el) > 1.0_pReal) then + call random_number(rnd) + sourceProbability(s,ipc,ip,el) = rnd + !$OMP FLUSH(sourceProbability) + endif + if (sourceProbability(s,ipc,ip,el) > 1.0_pReal - nSources(s)) then + rhoDotMultiplication(s,1:4) = sum(rhoSglOriginal(s,1:4) * abs(v(s,1:4))) / meshlength + endif + else + sourceProbability(s,ipc,ip,el) = 2.0_pReal + rhoDotMultiplication(s,1:4) = & + (sum(abs(gdot(s,1:2))) * fEdgeMultiplication(instance) + sum(abs(gdot(s,3:4)))) & + / burgers(s,instance) * sqrt(rhoForest(s)) / lambda0(s,instance) + endif + enddo +#ifndef _OPENMP + if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt & + .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)& + .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then + write(6,'(a,/,4(12x,12(f12.5,1x),/))') '<< CONST >> sources', nSources + write(6,*) + endif +#endif + else + rhoDotMultiplication(1:ns,1:4) = spread( & + (sum(abs(gdot(1:ns,1:2)),2) * fEdgeMultiplication(instance) + sum(abs(gdot(1:ns,3:4)),2)) & + * sqrt(rhoForest(1:ns)) / lambda0(1:ns,instance) / burgers(1:ns,instance), 2, 4) + endif +endif + + + +!**************************************************************************** +!*** calculate dislocation fluxes (only for nonlocal plasticity) + +rhoDotFlux = 0.0_pReal + +if (.not. phase_localPlasticity(material_phase(ipc,ip,el))) then ! only for nonlocal plasticity + + + !*** check CFL (Courant-Friedrichs-Lewy) condition for flux + + if (any( abs(gdot) > 0.0_pReal & ! any active slip system ... + .and. CFLfactor(instance) * abs(v) * timestep & + > mesh_ipVolume(ip,el) / maxval(mesh_ipArea(:,ip,el)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here) +#ifndef _OPENMP + if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt) then + write(6,'(a,i5,a,i2)') '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip + write(6,'(a,e10.3,a,e10.3)') '<< CONST >> velocity is at ', & + maxval(abs(v), abs(gdot) > 0.0_pReal & + .and. CFLfactor(instance) * abs(v) * timestep & + > mesh_ipVolume(ip,el) / maxval(mesh_ipArea(:,ip,el))), & + ' at a timestep of ',timestep + write(6,'(a)') '<< CONST >> enforcing cutback !!!' + endif +#endif + constitutive_nonlocal_dotState = DAMASK_NaN ! -> return NaN and, hence, enforce cutback + return + endif + + + !*** be aware of the definition of lattice_st = lattice_sd x lattice_sn !!! + !*** opposite sign to our p vector in the (s,p,n) triplet !!! + + m(1:3,1:ns,1) = lattice_sd(1:3, slipSystemLattice(1:ns,instance), phase) + m(1:3,1:ns,2) = -lattice_sd(1:3, slipSystemLattice(1:ns,instance), phase) + m(1:3,1:ns,3) = -lattice_st(1:3, slipSystemLattice(1:ns,instance), phase) + m(1:3,1:ns,4) = lattice_st(1:3, slipSystemLattice(1:ns,instance), phase) + + my_Fe = Fe(1:3,1:3,ipc,ip,el) + my_F = math_mul33x33(my_Fe, Fp(1:3,1:3,ipc,ip,el)) + + do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el)))) ! loop through my neighbors + neighbor_el = mesh_ipNeighborhood(1,n,ip,el) + neighbor_ip = mesh_ipNeighborhood(2,n,ip,el) + neighbor_n = mesh_ipNeighborhood(3,n,ip,el) + + opposite_neighbor = n + mod(n,2_pInt) - mod(n+1_pInt,2_pInt) + opposite_el = mesh_ipNeighborhood(1,opposite_neighbor,ip,el) + opposite_ip = mesh_ipNeighborhood(2,opposite_neighbor,ip,el) + opposite_n = mesh_ipNeighborhood(3,opposite_neighbor,ip,el) + + if (neighbor_n > 0_pInt) then ! if neighbor exists, average deformation gradient + neighbor_instance = phase_plasticityInstance(material_phase(ipc,neighbor_ip,neighbor_el)) + neighbor_Fe = Fe(1:3,1:3,ipc,neighbor_ip,neighbor_el) + neighbor_F = math_mul33x33(neighbor_Fe, Fp(1:3,1:3,ipc,neighbor_ip,neighbor_el)) + Favg = 0.5_pReal * (my_F + neighbor_F) + else ! if no neighbor, take my value as average + Favg = my_F + endif + + + !* FLUX FROM MY NEIGHBOR TO ME + !* This is only considered, if I have a neighbor of nonlocal plasticity (also nonlocal constitutive law with local properties) that is at least a little bit compatible. + !* If it's not at all compatible, no flux is arriving, because everything is dammed in front of my neighbor's interface. + !* The entering flux from my neighbor will be distributed on my slip systems according to the compatibility + + considerEnteringFlux = .false. + neighbor_v = 0.0_pReal ! needed for check of sign change in flux density below + neighbor_rhoSgl = 0.0_pReal + if (neighbor_n > 0_pInt) then + if (phase_plasticity(material_phase(1,neighbor_ip,neighbor_el)) == PLASTICITY_NONLOCAL_ID & + .and. any(compatibility(:,:,:,n,ip,el) > 0.0_pReal)) & + considerEnteringFlux = .true. + endif + + if (considerEnteringFlux) then + if(numerics_timeSyncing .and. (subfrac(ipc,neighbor_ip,neighbor_el) /= subfrac(ipc,ip,el))) then ! for timesyncing: in case of a timestep at the interface we have to use "state0" to make sure that fluxes n both sides are equal + forall (s = 1:ns, t = 1_pInt:4_pInt) + neighbor_v(s,t) = plasticState(mapping(2,ipc,neighbor_ip,neighbor_el))% & + state0(iV(s,t,neighbor_instance),mapping(1,ipc,ip,el)) + neighbor_rhoSgl(s,t) = max(plasticState(mapping(2,ipc,neighbor_ip,neighbor_el))% & + state0(iV(s,t,neighbor_instance),mapping(1,ipc,ip,el)), 0.0_pReal) + endforall + else + forall (s = 1:ns, t = 1_pInt:4_pInt) + neighbor_v(s,t) = plasticState(mapping(2,ipc,neighbor_ip,neighbor_el))% & + state(iV(s,t,neighbor_instance),mapping(1,ipc,ip,el)) + neighbor_rhoSgl(s,t) = max(plasticState(mapping(2,ipc,neighbor_ip,neighbor_el))% & + state(iV(s,t,neighbor_instance),mapping(1,ipc,ip,el)), 0.0_pReal) + endforall + endif + where (neighbor_rhoSgl * mesh_ipVolume(neighbor_ip,neighbor_el) ** 0.667_pReal < significantN(instance) & + .or. neighbor_rhoSgl < significantRho(instance)) & + neighbor_rhoSgl = 0.0_pReal + normal_neighbor2me_defConf = math_det33(Favg) * math_mul33x3(math_inv33(transpose(Favg)), & + mesh_ipAreaNormal(1:3,neighbor_n,neighbor_ip,neighbor_el)) ! calculate the normal of the interface in (average) deformed configuration (now pointing from my neighbor to me!!!) + normal_neighbor2me = math_mul33x3(transpose(neighbor_Fe), normal_neighbor2me_defConf) & + / math_det33(neighbor_Fe) ! interface normal in the lattice configuration of my neighbor + area = mesh_ipArea(neighbor_n,neighbor_ip,neighbor_el) * math_norm3(normal_neighbor2me) + normal_neighbor2me = normal_neighbor2me / math_norm3(normal_neighbor2me) ! normalize the surface normal to unit length + do s = 1_pInt,ns + do t = 1_pInt,4_pInt + c = (t + 1_pInt) / 2 + topp = t + mod(t,2_pInt) - mod(t+1_pInt,2_pInt) + if (neighbor_v(s,t) * math_mul3x3(m(1:3,s,t), normal_neighbor2me) > 0.0_pReal & ! flux from my neighbor to me == entering flux for me + .and. v(s,t) * neighbor_v(s,t) >= 0.0_pReal ) then ! ... only if no sign change in flux density + lineLength = neighbor_rhoSgl(s,t) * neighbor_v(s,t) & + * math_mul3x3(m(1:3,s,t), normal_neighbor2me) * area ! positive line length that wants to enter through this interface + where (compatibility(c,1_pInt:ns,s,n,ip,el) > 0.0_pReal) & ! positive compatibility... + rhoDotFlux(1_pInt:ns,t) = rhoDotFlux(1_pInt:ns,t) & + + lineLength / mesh_ipVolume(ip,el) & ! ... transferring to equally signed mobile dislocation type + * compatibility(c,1_pInt:ns,s,n,ip,el) ** 2.0_pReal + where (compatibility(c,1_pInt:ns,s,n,ip,el) < 0.0_pReal) & ! ..negative compatibility... + rhoDotFlux(1_pInt:ns,topp) = rhoDotFlux(1_pInt:ns,topp) & + + lineLength / mesh_ipVolume(ip,el) & ! ... transferring to opposite signed mobile dislocation type + * compatibility(c,1_pInt:ns,s,n,ip,el) ** 2.0_pReal + endif + enddo + enddo + endif + + + !* FLUX FROM ME TO MY NEIGHBOR + !* This is not considered, if my opposite neighbor has a different constitutive law than nonlocal (still considered for nonlocal law with lcal properties). + !* Then, we assume, that the opposite(!) neighbor sends an equal amount of dislocations to me. + !* So the net flux in the direction of my neighbor is equal to zero: + !* leaving flux to neighbor == entering flux from opposite neighbor + !* In case of reduced transmissivity, part of the leaving flux is stored as dead dislocation density. + !* That means for an interface of zero transmissivity the leaving flux is fully converted to dead dislocations. + + considerLeavingFlux = .true. + if (opposite_n > 0_pInt) then + if (phase_plasticity(material_phase(1,opposite_ip,opposite_el)) /= PLASTICITY_NONLOCAL_ID) & + considerLeavingFlux = .false. + endif + + if (considerLeavingFlux) then + + !* timeSyncing mode: If the central ip has zero subfraction, always use "state0". This is needed in case of + !* a synchronization step for the central ip, because then "state" contains the values at the end of the + !* previously converged full time step. Also, if either me or my neighbor has zero subfraction, we have to + !* use "state0" to make sure that fluxes on both sides of the (potential) timestep are equal. + my_rhoSgl = rhoSgl + my_v = v + if(numerics_timeSyncing) then + if (subfrac(ipc,ip,el) == 0.0_pReal) then + my_rhoSgl = rhoSgl0 + my_v = v0 + elseif (neighbor_n > 0_pInt) then + if (subfrac(ipc,neighbor_ip,neighbor_el) == 0.0_pReal) then + my_rhoSgl = rhoSgl0 + my_v = v0 + endif + endif + endif + + normal_me2neighbor_defConf = math_det33(Favg) & + * math_mul33x3(math_inv33(math_transpose33(Favg)), & + mesh_ipAreaNormal(1:3,n,ip,el)) ! calculate the normal of the interface in (average) deformed configuration (pointing from me to my neighbor!!!) + normal_me2neighbor = math_mul33x3(math_transpose33(my_Fe), normal_me2neighbor_defConf) & + / math_det33(my_Fe) ! interface normal in my lattice configuration + area = mesh_ipArea(n,ip,el) * math_norm3(normal_me2neighbor) + normal_me2neighbor = normal_me2neighbor / math_norm3(normal_me2neighbor) ! normalize the surface normal to unit length + do s = 1_pInt,ns + do t = 1_pInt,4_pInt + c = (t + 1_pInt) / 2_pInt + if (my_v(s,t) * math_mul3x3(m(1:3,s,t), normal_me2neighbor) > 0.0_pReal ) then ! flux from me to my neighbor == leaving flux for me (might also be a pure flux from my mobile density to dead density if interface not at all transmissive) + if (my_v(s,t) * neighbor_v(s,t) >= 0.0_pReal) then ! no sign change in flux density + transmissivity = sum(compatibility(c,1_pInt:ns,s,n,ip,el)**2.0_pReal) ! overall transmissivity from this slip system to my neighbor + else ! sign change in flux density means sign change in stress which does not allow for dislocations to arive at the neighbor + transmissivity = 0.0_pReal + endif + lineLength = my_rhoSgl(s,t) * my_v(s,t) & + * math_mul3x3(m(1:3,s,t), normal_me2neighbor) * area ! positive line length of mobiles that wants to leave through this interface + rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / mesh_ipVolume(ip,el) ! subtract dislocation flux from current type + rhoDotFlux(s,t+4_pInt) = rhoDotFlux(s,t+4_pInt) & + + lineLength / mesh_ipVolume(ip,el) * (1.0_pReal - transmissivity) & + * sign(1.0_pReal, my_v(s,t)) ! dislocation flux that is not able to leave through interface (because of low transmissivity) will remain as immobile single density at the material point + endif + enddo + enddo + endif + + enddo ! neighbor loop +endif + + + +!**************************************************************************** +!*** calculate dipole formation and annihilation + +!*** formation by glide + +do c = 1_pInt,2_pInt + rhoDotSingle2DipoleGlide(1:ns,2*c-1) = -2.0_pReal * dUpper(1:ns,c) / burgers(1:ns,instance) & + * (rhoSgl(1:ns,2*c-1) * abs(gdot(1:ns,2*c)) & ! negative mobile --> positive mobile + + rhoSgl(1:ns,2*c) * abs(gdot(1:ns,2*c-1)) & ! positive mobile --> negative mobile + + abs(rhoSgl(1:ns,2*c+4)) * abs(gdot(1:ns,2*c-1))) ! positive mobile --> negative immobile + + rhoDotSingle2DipoleGlide(1:ns,2*c) = -2.0_pReal * dUpper(1:ns,c) / burgers(1:ns,instance) & + * (rhoSgl(1:ns,2*c-1) * abs(gdot(1:ns,2*c)) & ! negative mobile --> positive mobile + + rhoSgl(1:ns,2*c) * abs(gdot(1:ns,2*c-1)) & ! positive mobile --> negative mobile + + abs(rhoSgl(1:ns,2*c+3)) * abs(gdot(1:ns,2*c))) ! negative mobile --> positive immobile + + rhoDotSingle2DipoleGlide(1:ns,2*c+3) = -2.0_pReal * dUpper(1:ns,c) / burgers(1:ns,instance) & + * rhoSgl(1:ns,2*c+3) * abs(gdot(1:ns,2*c)) ! negative mobile --> positive immobile + + rhoDotSingle2DipoleGlide(1:ns,2*c+4) = -2.0_pReal * dUpper(1:ns,c) / burgers(1:ns,instance) & + * rhoSgl(1:ns,2*c+4) * abs(gdot(1:ns,2*c-1)) ! positive mobile --> negative immobile + + rhoDotSingle2DipoleGlide(1:ns,c+8) = - rhoDotSingle2DipoleGlide(1:ns,2*c-1) & + - rhoDotSingle2DipoleGlide(1:ns,2*c) & + + abs(rhoDotSingle2DipoleGlide(1:ns,2*c+3)) & + + abs(rhoDotSingle2DipoleGlide(1:ns,2*c+4)) +enddo + + +!*** athermal annihilation + +rhoDotAthermalAnnihilation = 0.0_pReal + +forall (c=1_pInt:2_pInt) & + rhoDotAthermalAnnihilation(1:ns,c+8_pInt) = -2.0_pReal * dLower(1:ns,c) / burgers(1:ns,instance) & + * ( 2.0_pReal * (rhoSgl(1:ns,2*c-1) * abs(gdot(1:ns,2*c)) + rhoSgl(1:ns,2*c) * abs(gdot(1:ns,2*c-1))) & ! was single hitting single + + 2.0_pReal * (abs(rhoSgl(1:ns,2*c+3)) * abs(gdot(1:ns,2*c)) + abs(rhoSgl(1:ns,2*c+4)) * abs(gdot(1:ns,2*c-1))) & ! was single hitting immobile single or was immobile single hit by single + + rhoDip(1:ns,c) * (abs(gdot(1:ns,2*c-1)) + abs(gdot(1:ns,2*c)))) ! single knocks dipole constituent +! annihilated screw dipoles leave edge jogs behind on the colinear system +if (lattice_structure(phase) == LATTICE_fcc_ID) then ! only fcc + forall (s = 1:ns, colinearSystem(s,instance) > 0_pInt) & + rhoDotAthermalAnnihilation(colinearSystem(s,instance),1:2) = - rhoDotAthermalAnnihilation(s,10) & + * 0.25_pReal * sqrt(rhoForest(s)) * (dUpper(s,2) + dLower(s,2)) * edgeJogFactor(instance) +endif + + +!*** thermally activated annihilation of edge dipoles by climb + +rhoDotThermalAnnihilation = 0.0_pReal +selfDiffusion = Dsd0(instance) * exp(-selfDiffusionEnergy(instance) / (KB * Temperature)) +vClimb = atomicVolume(instance) * selfDiffusion / ( KB * Temperature ) & + * lattice_mu(phase) / ( 2.0_pReal * PI * (1.0_pReal-lattice_nu(phase)) ) & + * 2.0_pReal / ( dUpper(1:ns,1) + dLower(1:ns,1) ) +forall (s = 1_pInt:ns, dUpper(s,1) > dLower(s,1)) & + rhoDotThermalAnnihilation(s,9) = max(- 4.0_pReal * rhoDip(s,1) * vClimb(s) / (dUpper(s,1) - dLower(s,1)), & + - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) & + - rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have + + + +!**************************************************************************** +!*** assign the rates of dislocation densities to my dotState +!*** if evolution rates lead to negative densities, a cutback is enforced + +rhoDot = 0.0_pReal +rhoDot = rhoDotFlux & + + rhoDotMultiplication & + + rhoDotSingle2DipoleGlide & + + rhoDotAthermalAnnihilation & + + rhoDotThermalAnnihilation + +if (numerics_integrationMode == 1_pInt) then ! save rates for output if in central integration mode + rhoDotFluxOutput(1:ns,1:8,ipc,ip,el) = rhoDotFlux(1:ns,1:8) + rhoDotMultiplicationOutput(1:ns,1:2,ipc,ip,el) = rhoDotMultiplication(1:ns,[1,3]) + rhoDotSingle2DipoleGlideOutput(1:ns,1:2,ipc,ip,el) = rhoDotSingle2DipoleGlide(1:ns,9:10) + rhoDotAthermalAnnihilationOutput(1:ns,1:2,ipc,ip,el) = rhoDotAthermalAnnihilation(1:ns,9:10) + rhoDotThermalAnnihilationOutput(1:ns,1:2,ipc,ip,el) = rhoDotThermalAnnihilation(1:ns,9:10) + rhoDotEdgeJogsOutput(1:ns,ipc,ip,el) = 2.0_pReal * rhoDotThermalAnnihilation(1:ns,1) +endif + + +#ifndef _OPENMP + if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt & + .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)& + .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then + write(6,'(a,/,4(12x,12(e12.5,1x),/))') '<< CONST >> dislocation multiplication', & + rhoDotMultiplication(1:ns,1:4) * timestep + write(6,'(a,/,8(12x,12(e12.5,1x),/))') '<< CONST >> dislocation flux', & + rhoDotFlux(1:ns,1:8) * timestep + write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> dipole formation by glide', & + rhoDotSingle2DipoleGlide * timestep + write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> athermal dipole annihilation', & + rhoDotAthermalAnnihilation * timestep + write(6,'(a,/,2(12x,12(e12.5,1x),/))') '<< CONST >> thermally activated dipole annihilation', & + rhoDotThermalAnnihilation(1:ns,9:10) * timestep + write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> total density change', & + rhoDot * timestep + write(6,'(a,/,10(12x,12(f12.5,1x),/))') '<< CONST >> relative density change', & + rhoDot(1:ns,1:8) * timestep / (abs(rhoSglOriginal)+1.0e-10), & + rhoDot(1:ns,9:10) * timestep / (rhoDipOriginal+1.0e-10) + write(6,*) + endif +#endif + + +if ( any(rhoSglOriginal(1:ns,1:4) + rhoDot(1:ns,1:4) * timestep < -aTolRho(instance)) & + .or. any(rhoDipOriginal(1:ns,1:2) + rhoDot(1:ns,9:10) * timestep < -aTolRho(instance))) then +#ifndef _OPENMP + if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt) then + write(6,'(a,i5,a,i2)') '<< CONST >> evolution rate leads to negative density at el ',el,' ip ',ip + write(6,'(a)') '<< CONST >> enforcing cutback !!!' + endif +#endif + constitutive_nonlocal_dotState = DAMASK_NaN + return +else + forall (s = 1:ns, t = 1_pInt:4_pInt) + constitutive_nonlocal_dotState(iRhoU(s,t,instance)) = rhoDot(s,t) + constitutive_nonlocal_dotState(iRhoB(s,t,instance)) = rhoDot(s,t+4_pInt) + endforall + forall (s = 1:ns, c = 1_pInt:2_pInt) & + constitutive_nonlocal_dotState(iRhoD(s,c,instance)) = rhoDot(s,c+8_pInt) + forall (s = 1:ns) & + constitutive_nonlocal_dotState(iGamma(s,instance)) = sum(gdot(s,1:4)) +endif + +end function constitutive_nonlocal_dotState + + + +#else !-------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !-------------------------------------------------------------------------------------------------- @@ -2759,7 +4311,7 @@ endif end function constitutive_nonlocal_dotState - +#endif !********************************************************************* !* COMPATIBILITY UPDATE * @@ -2932,8 +4484,365 @@ enddo ! neighbor cycle compatibility(1:2,1:ns,1:ns,1:Nneighbors,i,e) = my_compatibility end subroutine constitutive_nonlocal_updateCompatibility +#ifdef NEWSTATE +!********************************************************************* +!* calculates quantities characterizing the microstructure * +!********************************************************************* +function constitutive_nonlocal_dislocationstress(mapp, Fe, ipc, ip, el) + +use math, only: math_mul33x33, & + math_mul33x3, & + math_invert33, & + math_transpose33, & + pi +use mesh, only: mesh_NcpElems, & + mesh_maxNips, & + mesh_element, & + mesh_node0, & + mesh_cellCenterCoordinates, & + mesh_ipVolume, & + mesh_periodicSurface, & + FE_Nips, & + FE_geomtype +use material, only: homogenization_maxNgrains, & + material_phase, & + phase_localPlasticity, & + phase_plasticityInstance, & + plasticState +use lattice, only: lattice_mu, & + lattice_nu + +implicit none + +!*** input variables +integer(pInt), intent(in) :: ipc, & !< current grain ID + ip, & !< current integration point + el !< current element +real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + Fe !< elastic deformation gradient +!type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & +! state !< microstructural state +integer(pInt),dimension(2,homogenization_maxngrains,mesh_maxNips,mesh_ncpelems), intent(in) :: & + mapp +!*** output variables +real(pReal), dimension(3,3) :: constitutive_nonlocal_dislocationstress + +!*** local variables +integer(pInt) neighbor_el, & !< element number of neighbor material point + neighbor_ip, & !< integration point of neighbor material point + instance, & !< my instance of this plasticity + neighbor_instance, & !< instance of this plasticity of neighbor material point + phase, & + neighbor_phase, & + ns, & !< total number of active slip systems at my material point + neighbor_ns, & !< total number of active slip systems at neighbor material point + c, & !< index of dilsocation character (edge, screw) + s, & !< slip system index + t, & !< index of dilsocation type (e+, e-, s+, s-, used e+, used e-, used s+, used s-) + dir, & + deltaX, deltaY, deltaZ, & + side, & + j +integer(pInt), dimension(2,3) :: periodicImages +real(pReal) x, y, z, & !< coordinates of connection vector in neighbor lattice frame + xsquare, ysquare, zsquare, & !< squares of respective coordinates + distance, & !< length of connection vector + segmentLength, & !< segment length of dislocations + lambda, & + R, Rsquare, Rcube, & + denominator, & + flipSign, & + neighbor_ipVolumeSideLength, & + detFe +real(pReal), dimension(3) :: connection, & !< connection vector between me and my neighbor in the deformed configuration + connection_neighborLattice, & !< connection vector between me and my neighbor in the lattice configuration of my neighbor + connection_neighborSlip, & !< connection vector between me and my neighbor in the slip system frame of my neighbor + maxCoord, minCoord, & + meshSize, & + coords, & !< x,y,z coordinates of cell center of ip volume + neighbor_coords !< x,y,z coordinates of cell center of neighbor ip volume +real(pReal), dimension(3,3) :: sigma, & !< dislocation stress for one slip system in neighbor material point's slip system frame + Tdislo_neighborLattice, & !< dislocation stress as 2nd Piola-Kirchhoff stress at neighbor material point + invFe, & !< inverse of my elastic deformation gradient + neighbor_invFe, & + neighborLattice2myLattice !< mapping from neighbor MPs lattice configuration to my lattice configuration +real(pReal), dimension(2,2,maxval(totalNslip)) :: & + neighbor_rhoExcess !< excess density at neighbor material point (edge/screw,mobile/dead,slipsystem) +real(pReal), dimension(2,maxval(totalNslip)) :: & + rhoExcessDead +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),8) :: & + rhoSgl ! single dislocation density (edge+, edge-, screw+, screw-, used edge+, used edge-, used screw+, used screw-) +logical inversionError + +phase = material_phase(ipc,ip,el) +instance = phase_plasticityInstance(phase) +ns = totalNslip(instance) + +!*** get basic states + +forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) +! rhoSgl(s,t) = max(state(ipc,ip,el)%p(iRhoU(s,t,instance)), 0.0_pReal) ! ensure positive single mobile densities +! rhoSgl(s,t+4_pInt) = state(ipc,ip,el)%p(iRhoB(s,t,instance)) + + rhoSgl(s,t) = max(plasticState(mapp(2,ipc,ip,el))%state(iRhoU(s,t,instance),mapp(1,ipc,ip,el)), & + 0.0_pReal) ! ensure positive single mobile densities + rhoSgl(s,t+4_pInt) = plasticState(mapp(2,ipc,ip,el))%state(iRhoB(s,t,instance),mapp(1,ipc,ip,el)) +endforall + +!*** calculate the dislocation stress of the neighboring excess dislocation densities +!*** zero for material points of local plasticity + +constitutive_nonlocal_dislocationstress = 0.0_pReal + +if (.not. phase_localPlasticity(phase)) then + call math_invert33(Fe(1:3,1:3,ipc,ip,el), invFe, detFe, inversionError) + + !* in case of periodic surfaces we have to find out how many periodic images in each direction we need + + do dir = 1_pInt,3_pInt + maxCoord(dir) = maxval(mesh_node0(dir,:)) + minCoord(dir) = minval(mesh_node0(dir,:)) + enddo + meshSize = maxCoord - minCoord + coords = mesh_cellCenterCoordinates(ip,el) + periodicImages = 0_pInt + do dir = 1_pInt,3_pInt + if (mesh_periodicSurface(dir)) then + periodicImages(1,dir) = floor((coords(dir) - cutoffRadius(instance) - minCoord(dir)) / meshSize(dir), pInt) + periodicImages(2,dir) = ceiling((coords(dir) + cutoffRadius(instance) - maxCoord(dir)) / meshSize(dir), pInt) + endif + enddo + + + !* loop through all material points (also through their periodic images if present), + !* but only consider nonlocal neighbors within a certain cutoff radius R + + do neighbor_el = 1_pInt,mesh_NcpElems +ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el))) + neighbor_phase = material_phase(ipc,neighbor_ip,neighbor_el) + if (phase_localPlasticity(neighbor_phase)) then + cycle + endif + neighbor_instance = phase_plasticityInstance(neighbor_phase) + neighbor_ns = totalNslip(neighbor_instance) + call math_invert33(Fe(1:3,1:3,1,neighbor_ip,neighbor_el), neighbor_invFe, detFe, inversionError) + neighbor_ipVolumeSideLength = mesh_ipVolume(neighbor_ip,neighbor_el) ** (1.0_pReal/3.0_pReal) ! reference volume used here + forall (s = 1_pInt:neighbor_ns, c = 1_pInt:2_pInt) + neighbor_rhoExcess(c,1,s) = plasticState(mapp(2,ipc,neighbor_ip,neighbor_el))% & + state(iRhoU(s,2*c-1,neighbor_instance),mapp(1,ipc,neighbor_ip,neighbor_el)) & ! positive mobiles + - plasticState(mapp(2,ipc,neighbor_ip,neighbor_el))% & + state(iRhoU(s,2*c,neighbor_instance),mapp(1,ipc,neighbor_ip,neighbor_el)) ! negative mobiles + neighbor_rhoExcess(c,2,s) = abs(plasticState(mapp(2,ipc,neighbor_ip,neighbor_el))% & + state(iRhoB(s,2*c-1,neighbor_instance),mapp(1,ipc,neighbor_ip,neighbor_el))) & ! positive deads + - abs(plasticState(mapp(2,ipc,neighbor_ip,neighbor_el))% & + state(iRhoB(s,2*c,neighbor_instance),mapp(1,ipc,neighbor_ip,neighbor_el))) ! negative deads + endforall + Tdislo_neighborLattice = 0.0_pReal + do deltaX = periodicImages(1,1),periodicImages(2,1) + do deltaY = periodicImages(1,2),periodicImages(2,2) + do deltaZ = periodicImages(1,3),periodicImages(2,3) + + + !* regular case + + if (neighbor_el /= el .or. neighbor_ip /= ip & + .or. deltaX /= 0_pInt .or. deltaY /= 0_pInt .or. deltaZ /= 0_pInt) then + + neighbor_coords = mesh_cellCenterCoordinates(neighbor_ip,neighbor_el) & + + [real(deltaX,pReal), real(deltaY,pReal), real(deltaZ,pReal)] * meshSize + connection = neighbor_coords - coords + distance = sqrt(sum(connection * connection)) + if (distance > cutoffRadius(instance)) then + cycle + endif + + + !* the segment length is the minimum of the third root of the control volume and the ip distance + !* this ensures, that the central MP never sits on a neighbor dislocation segment + + connection_neighborLattice = math_mul33x3(neighbor_invFe, connection) + segmentLength = min(neighbor_ipVolumeSideLength, distance) + + + !* loop through all slip systems of the neighbor material point + !* and add up the stress contributions from egde and screw excess on these slip systems (if significant) + + do s = 1_pInt,neighbor_ns + if (all(abs(neighbor_rhoExcess(:,:,s)) < significantRho(instance))) cycle ! not significant + + + !* map the connection vector from the lattice into the slip system frame + + connection_neighborSlip = math_mul33x3(lattice2slip(1:3,1:3,s,neighbor_instance), & + connection_neighborLattice) + + + !* edge contribution to stress + sigma = 0.0_pReal + + x = connection_neighborSlip(1) + y = connection_neighborSlip(2) + z = connection_neighborSlip(3) + xsquare = x * x + ysquare = y * y + zsquare = z * z + + do j = 1_pInt,2_pInt + if (abs(neighbor_rhoExcess(1,j,s)) < significantRho(instance)) then + cycle + elseif (j > 1_pInt) then +! x = connection_neighborSlip(1) & +! + sign(0.5_pReal * segmentLength, & +! state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,1,neighbor_instance)) & +! - state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,2,neighbor_instance))) + x = connection_neighborSlip(1) & + + sign(0.5_pReal * segmentLength, & + plasticState(mapp(2,ipc,neighbor_ip,neighbor_el))% & + state(iRhoB(s,1,neighbor_instance),mapp(1,ipc,neighbor_ip,neighbor_el)) & + - plasticState(mapp(2,ipc,neighbor_ip,neighbor_el))% & + state(iRhoB(s,2,neighbor_instance),mapp(1,ipc,neighbor_ip,neighbor_el))) + xsquare = x * x + endif + + flipSign = sign(1.0_pReal, -y) + do side = 1_pInt,-1_pInt,-2_pInt + lambda = real(side,pReal) * 0.5_pReal * segmentLength - y + R = sqrt(xsquare + zsquare + lambda * lambda) + Rsquare = R * R + Rcube = Rsquare * R + denominator = R * (R + flipSign * lambda) + if (denominator == 0.0_pReal) exit ipLoop + + sigma(1,1) = sigma(1,1) - real(side,pReal) & + * flipSign * z / denominator & + * (1.0_pReal + xsquare / Rsquare + xsquare / denominator) & + * neighbor_rhoExcess(1,j,s) + sigma(2,2) = sigma(2,2) - real(side,pReal) & + * (flipSign * 2.0_pReal * lattice_nu(phase) * z / denominator + z * lambda / Rcube) & + * neighbor_rhoExcess(1,j,s) + sigma(3,3) = sigma(3,3) + real(side,pReal) & + * flipSign * z / denominator & + * (1.0_pReal - zsquare / Rsquare - zsquare / denominator) & + * neighbor_rhoExcess(1,j,s) + sigma(1,2) = sigma(1,2) + real(side,pReal) & + * x * z / Rcube * neighbor_rhoExcess(1,j,s) + sigma(1,3) = sigma(1,3) + real(side,pReal) & + * flipSign * x / denominator & + * (1.0_pReal - zsquare / Rsquare - zsquare / denominator) & + * neighbor_rhoExcess(1,j,s) + sigma(2,3) = sigma(2,3) - real(side,pReal) & + * (lattice_nu(phase) / R - zsquare / Rcube) * neighbor_rhoExcess(1,j,s) + enddo + enddo + + !* screw contribution to stress + + x = connection_neighborSlip(1) ! have to restore this value, because position might have been adapted for edge deads before + do j = 1_pInt,2_pInt + if (abs(neighbor_rhoExcess(2,j,s)) < significantRho(instance)) then + cycle + elseif (j > 1_pInt) then + y = connection_neighborSlip(2) & + + sign(0.5_pReal * segmentLength, & + plasticState(mapp(2,ipc,neighbor_ip,neighbor_el))% & + state(iRhoB(s,3,neighbor_instance),mapp(1,ipc,neighbor_ip,neighbor_el)) & + - plasticState(mapp(2,ipc,neighbor_ip,neighbor_el))% & + state(iRhoB(s,4,neighbor_instance),mapp(1,ipc,neighbor_ip,neighbor_el))) + ysquare = y * y + endif + + flipSign = sign(1.0_pReal, x) + do side = 1_pInt,-1_pInt,-2_pInt + lambda = x + real(side,pReal) * 0.5_pReal * segmentLength + R = sqrt(ysquare + zsquare + lambda * lambda) + Rsquare = R * R + Rcube = Rsquare * R + denominator = R * (R + flipSign * lambda) + if (denominator == 0.0_pReal) exit ipLoop + + sigma(1,2) = sigma(1,2) - real(side,pReal) * flipSign * z & + * (1.0_pReal - lattice_nu(phase)) / denominator & + * neighbor_rhoExcess(2,j,s) + sigma(1,3) = sigma(1,3) + real(side,pReal) * flipSign * y & + * (1.0_pReal - lattice_nu(phase)) / denominator & + * neighbor_rhoExcess(2,j,s) + enddo + enddo + + if (all(abs(sigma) < 1.0e-10_pReal)) cycle ! SIGMA IS NOT A REAL STRESS, THATS WHY WE NEED A REALLY SMALL VALUE HERE + + !* copy symmetric parts + + sigma(2,1) = sigma(1,2) + sigma(3,1) = sigma(1,3) + sigma(3,2) = sigma(2,3) + + + !* scale stresses and map them into the neighbor material point's lattice configuration + + sigma = sigma * lattice_mu(neighbor_phase) * burgers(s,neighbor_instance) & + / (4.0_pReal * pi * (1.0_pReal - lattice_nu(neighbor_phase))) & + * mesh_ipVolume(neighbor_ip,neighbor_el) / segmentLength ! reference volume is used here (according to the segment length calculation) + Tdislo_neighborLattice = Tdislo_neighborLattice & + + math_mul33x33(math_transpose33(lattice2slip(1:3,1:3,s,neighbor_instance)), & + math_mul33x33(sigma, lattice2slip(1:3,1:3,s,neighbor_instance))) + + enddo ! slip system loop + + + !* special case of central ip volume + !* only consider dead dislocations + !* we assume that they all sit at a distance equal to half the third root of V + !* in direction of the according slip direction + + else + + forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) & + rhoExcessDead(c,s) = plasticState(mapp(2,ipc,ip,el))% & + state(iRhoB(s,2*c-1,instance),mapp(1,ipc,ip,el)) & ! positive deads (here we use symmetry: if this has negative sign it is treated as negative density at positive position instead of positive density at negative position) + + plasticState(mapp(2,ipc,neighbor_ip,neighbor_el))% & + state(iRhoB(s,2*c,instance),mapp(1,ipc,ip,el)) ! negative deads (here we use symmetry: if this has negative sign it is treated as positive density at positive position instead of negative density at negative position) + + do s = 1_pInt,ns + if (all(abs(rhoExcessDead(:,s)) < significantRho(instance))) cycle ! not significant + sigma = 0.0_pReal ! all components except for sigma13 are zero + sigma(1,3) = - (rhoExcessDead(1,s) + rhoExcessDead(2,s) * (1.0_pReal - lattice_nu(phase))) & + * neighbor_ipVolumeSideLength * lattice_mu(phase) * burgers(s,instance) & + / (sqrt(2.0_pReal) * pi * (1.0_pReal - lattice_nu(phase))) + sigma(3,1) = sigma(1,3) + + Tdislo_neighborLattice = Tdislo_neighborLattice & + + math_mul33x33(math_transpose33(lattice2slip(1:3,1:3,s,instance)), & + math_mul33x33(sigma, lattice2slip(1:3,1:3,s,instance))) + + enddo ! slip system loop + + endif + + enddo ! deltaZ loop + enddo ! deltaY loop + enddo ! deltaX loop + + + !* map the stress from the neighbor MP's lattice configuration into the deformed configuration + !* and back into my lattice configuration + + neighborLattice2myLattice = math_mul33x33(invFe, Fe(1:3,1:3,1,neighbor_ip,neighbor_el)) + constitutive_nonlocal_dislocationstress = constitutive_nonlocal_dislocationstress & + + math_mul33x33(neighborLattice2myLattice, & + math_mul33x33(Tdislo_neighborLattice, & + math_transpose33(neighborLattice2myLattice))) + + enddo ipLoop + enddo ! element loop + +endif + +end function constitutive_nonlocal_dislocationstress + +#else !********************************************************************* !* calculates quantities characterizing the microstructure * !********************************************************************* @@ -3272,8 +5181,526 @@ ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el)) endif end function constitutive_nonlocal_dislocationstress +#endif + +#ifdef NEWSTATE +!-------------------------------------------------------------------------------------------------- +!> @brief return array of constitutive results +!-------------------------------------------------------------------------------------------------- +function constitutive_nonlocal_postResults(Tstar_v,Fe,mapping,ipc,ip,el) + use math, only: & + math_mul6x6, & + math_mul33x3, & + math_mul33x33, & + pi + use mesh, only: & + mesh_NcpElems, & + mesh_maxNips + use material, only: & + homogenization_maxNgrains, & + material_phase, & + phase_plasticityInstance, & + phase_Noutput, & + plasticState + use lattice, only: & + lattice_Sslip_v, & + lattice_sd, & + lattice_st, & + lattice_sn, & + lattice_mu, & + lattice_nu + + implicit none + real(pReal), dimension(6), intent(in) :: & + Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation + real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + Fe !< elastic deformation gradient +! type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & +! state !< microstructure state +! type(p_vec), intent(in) :: & +! dotState ! evolution rate of microstructural state + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element +integer(pInt),dimension(2,homogenization_maxngrains,mesh_maxNips,mesh_ncpelems), intent(in) :: & + mapping + real(pReal), dimension(constitutive_nonlocal_sizePostResults(& + phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & + constitutive_nonlocal_postResults + + integer(pInt) :: & + phase, & + instance, & !< current instance of this plasticity + ns, & !< short notation for the total number of active slip systems + c, & !< character of dislocation + cs, & !< constitutive result index + o, & !< index of current output + t, & !< type of dislocation + s, & !< index of my current slip system + sLattice !< index of my current slip system according to lattice order + real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),8) :: & + rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles) + rhoDotSgl !< evolution rate of single dislocation densities (positive/negative screw and edge without dipoles) + real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),4) :: & + gdot, & !< shear rates + v !< velocities + real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & + rhoForest, & !< forest dislocation density + tauThreshold, & !< threshold shear stress + tau, & !< current resolved shear stress + tauBack !< back stress from pileups on same slip system + real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),2) :: & + rhoDip, & !< current dipole dislocation densities (screw and edge dipoles) + rhoDotDip, & !< evolution rate of dipole dislocation densities (screw and edge dipoles) + dLower, & !< minimum stable dipole distance for edges and screws + dUpper !< current maximum stable dipole distance for edges and screws + real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),2) :: & + m, & !< direction of dislocation motion for edge and screw (unit vector) + m_currentconf !< direction of dislocation motion for edge and screw (unit vector) in current configuration + real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & + n_currentconf !< slip system normal (unit vector) in current configuration + real(pReal), dimension(3,3) :: & + sigma + +phase = material_phase(ipc,ip,el) +instance = phase_plasticityInstance(phase) +ns = totalNslip(instance) + +cs = 0_pInt +constitutive_nonlocal_postResults = 0.0_pReal +!* short hand notations for state variables + +forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) + rhoSgl(s,t) = plasticState(mapping(2,ipc,ip,el))%state(iRhoU(s,t,instance),mapping(1,ipc,ip,el)) + rhoSgl(s,t+4_pInt) = plasticState(mapping(2,ipc,ip,el))%state(iRhoB(s,t,instance),mapping(1,ipc,ip,el)) + v(s,t) = plasticState(mapping(2,ipc,ip,el))%state(iV(s,t,instance),mapping(1,ipc,ip,el)) + rhoDotSgl(s,t) = plasticState(mapping(2,ipc,ip,el))%dotState(iRhoU(s,t,instance),mapping(1,ipc,ip,el)) + rhoDotSgl(s,t+4_pInt) = plasticState(mapping(2,ipc,ip,el))%dotState(iRhoB(s,t,instance),mapping(1,ipc,ip,el)) +endforall +forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) +! rhoDip(s,c) = state(ipc,ip,el)%p(iRhoD(s,c,instance)) +! rhoDotDip(s,c) = dotState%p(iRhoD(s,c,instance)) + + rhoDip(s,c) = plasticState(mapping(2,ipc,ip,el))%State(iRhoD(s,c,instance),mapping(1,ipc,ip,el)) + rhoDotDip(s,c) = plasticState(mapping(2,ipc,ip,el))%State(iRhoD(s,c,instance),mapping(1,ipc,ip,el)) +endforall +!tauBack = state(ipc,ip,el)%p(iTauB(1:ns,instance)) + +rhoForest = plasticState(mapping(2,ipc,ip,el))%state(iRhoF(1:ns,instance),mapping(1,ipc,ip,el)) +tauThreshold = plasticState(mapping(2,ipc,ip,el))%state(iTauF(1:ns,instance),mapping(1,ipc,ip,el)) +tauBack = plasticState(mapping(2,ipc,ip,el))%state(iTauB(1:ns,instance),mapping(1,ipc,ip,el)) + + + +!* Calculate shear rate + +forall (t = 1_pInt:4_pInt) & + gdot(1:ns,t) = rhoSgl(1:ns,t) * burgers(1:ns,instance) * v(1:ns,t) + + +!* calculate limits for stable dipole height + +do s = 1_pInt,ns + sLattice = slipSystemLattice(s,instance) + tau(s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,phase)) + tauBack(s) + if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal +enddo + +dLower = minDipoleHeight(1:ns,1:2,instance) +dUpper(1:ns,1) = lattice_mu(phase) * burgers(1:ns,instance) & + / (8.0_pReal * pi * (1.0_pReal - lattice_nu(phase)) * abs(tau)) +dUpper(1:ns,2) = lattice_mu(phase) * burgers(1:ns,instance) & + / (4.0_pReal * pi * abs(tau)) +forall (c = 1_pInt:2_pInt) & + dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) & + + abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), & + dUpper(1:ns,c)) +dUpper = max(dUpper,dLower) + +!*** dislocation motion + +m(1:3,1:ns,1) = lattice_sd(1:3,slipSystemLattice(1:ns,instance),phase) +m(1:3,1:ns,2) = -lattice_st(1:3,slipSystemLattice(1:ns,instance),phase) +forall (c = 1_pInt:2_pInt, s = 1_pInt:ns) & + m_currentconf(1:3,s,c) = math_mul33x3(Fe(1:3,1:3,ipc,ip,el), m(1:3,s,c)) +forall (s = 1_pInt:ns) & + n_currentconf(1:3,s) = math_mul33x3(Fe(1:3,1:3,ipc,ip,el), & + lattice_sn(1:3,slipSystemLattice(s,instance),phase)) + + +outputsLoop: do o = 1_pInt,phase_Noutput(material_phase(ipc,ip,el)) + select case(constitutive_nonlocal_outputID(o,instance)) + case (rho_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(abs(rhoSgl),2) + sum(rhoDip,2) + cs = cs + ns + + case (rho_sgl_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(abs(rhoSgl),2) + cs = cs + ns + + case (rho_sgl_mobile_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(abs(rhoSgl(1:ns,1:4)),2) + cs = cs + ns + + case (rho_sgl_immobile_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoSgl(1:ns,5:8),2) + cs = cs + ns + + case (rho_dip_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDip,2) + cs = cs + ns + + case (rho_edge_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(abs(rhoSgl(1:ns,[1,2,5,6])),2) + rhoDip(1:ns,1) + cs = cs + ns + + case (rho_sgl_edge_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(abs(rhoSgl(1:ns,[1,2,5,6])),2) + cs = cs + ns + + case (rho_sgl_edge_mobile_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoSgl(1:ns,1:2),2) + cs = cs + ns + + case (rho_sgl_edge_immobile_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoSgl(1:ns,5:6),2) + cs = cs + ns + + case (rho_sgl_edge_pos_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,1) + abs(rhoSgl(1:ns,5)) + cs = cs + ns + + case (rho_sgl_edge_pos_mobile_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,1) + cs = cs + ns + + case (rho_sgl_edge_pos_immobile_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,5) + cs = cs + ns + + case (rho_sgl_edge_neg_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,2) + abs(rhoSgl(1:ns,6)) + cs = cs + ns + + case (rho_sgl_edge_neg_mobile_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,2) + cs = cs + ns + + case (rho_sgl_edge_neg_immobile_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,6) + cs = cs + ns + + case (rho_dip_edge_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDip(1:ns,1) + cs = cs + ns + + case (rho_screw_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(abs(rhoSgl(1:ns,[3,4,7,8])),2) + rhoDip(1:ns,2) + cs = cs + ns + + case (rho_sgl_screw_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(abs(rhoSgl(1:ns,[3,4,7,8])),2) + cs = cs + ns + + case (rho_sgl_screw_mobile_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoSgl(1:ns,3:4),2) + cs = cs + ns + + case (rho_sgl_screw_immobile_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoSgl(1:ns,7:8),2) + cs = cs + ns + + case (rho_sgl_screw_pos_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,3) + abs(rhoSgl(1:ns,7)) + cs = cs + ns + + case (rho_sgl_screw_pos_mobile_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,3) + cs = cs + ns + + case (rho_sgl_screw_pos_immobile_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,7) + cs = cs + ns + + case (rho_sgl_screw_neg_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,4) + abs(rhoSgl(1:ns,8)) + cs = cs + ns + + case (rho_sgl_screw_neg_mobile_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,4) + cs = cs + ns + + case (rho_sgl_screw_neg_immobile_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,8) + cs = cs + ns + + case (rho_dip_screw_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDip(1:ns,2) + cs = cs + ns + + case (excess_rho_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = (rhoSgl(1:ns,1) + abs(rhoSgl(1:ns,5))) & + - (rhoSgl(1:ns,2) + abs(rhoSgl(1:ns,6))) & + + (rhoSgl(1:ns,3) + abs(rhoSgl(1:ns,7))) & + - (rhoSgl(1:ns,4) + abs(rhoSgl(1:ns,8))) + cs = cs + ns + + case (excess_rho_edge_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = (rhoSgl(1:ns,1) + abs(rhoSgl(1:ns,5))) & + - (rhoSgl(1:ns,2) + abs(rhoSgl(1:ns,6))) + cs = cs + ns + + case (excess_rho_screw_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = (rhoSgl(1:ns,3) + abs(rhoSgl(1:ns,7))) & + - (rhoSgl(1:ns,4) + abs(rhoSgl(1:ns,8))) + cs = cs + ns + + case (rho_forest_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoForest + cs = cs + ns + + case (delta_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = 1.0_pReal / sqrt(sum(abs(rhoSgl),2) + sum(rhoDip,2)) + cs = cs + ns + + case (delta_sgl_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = 1.0_pReal / sqrt(sum(abs(rhoSgl),2)) + cs = cs + ns + + case (delta_dip_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = 1.0_pReal / sqrt(sum(rhoDip,2)) + cs = cs + ns + + case (shearrate_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(gdot,2) + cs = cs + ns + + case (resolvedstress_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = tau + cs = cs + ns + + case (resolvedstress_back_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = tauBack + cs = cs + ns + + case (resolvedstress_external_ID) + do s = 1_pInt,ns + sLattice = slipSystemLattice(s,instance) + constitutive_nonlocal_postResults(cs+s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,phase)) + enddo + cs = cs + ns + + case (resistance_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = tauThreshold + cs = cs + ns + + case (rho_dot_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotSgl(1:ns,1:4),2) & + + sum(rhoDotSgl(1:ns,5:8)*sign(1.0_pReal,rhoSgl(1:ns,5:8)),2) & + + sum(rhoDotDip,2) + cs = cs + ns + + case (rho_dot_sgl_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotSgl(1:ns,1:4),2) & + + sum(rhoDotSgl(1:ns,5:8)*sign(1.0_pReal,rhoSgl(1:ns,5:8)),2) + cs = cs + ns + + case (rho_dot_sgl_mobile_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotSgl(1:ns,1:4),2) + cs = cs + ns + + case (rho_dot_dip_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotDip,2) + cs = cs + ns + + case (rho_dot_gen_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotMultiplicationOutput(1:ns,1,ipc,ip,el) & + + rhoDotMultiplicationOutput(1:ns,2,ipc,ip,el) + cs = cs + ns + + case (rho_dot_gen_edge_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotMultiplicationOutput(1:ns,1,ipc,ip,el) + cs = cs + ns + + case (rho_dot_gen_screw_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotMultiplicationOutput(1:ns,2,ipc,ip,el) + cs = cs + ns + + case (rho_dot_sgl2dip_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotSingle2DipoleGlideOutput(1:ns,1,ipc,ip,el) & + + rhoDotSingle2DipoleGlideOutput(1:ns,2,ipc,ip,el) + cs = cs + ns + + case (rho_dot_sgl2dip_edge_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotSingle2DipoleGlideOutput(1:ns,1,ipc,ip,el) + cs = cs + ns + + case (rho_dot_sgl2dip_screw_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotSingle2DipoleGlideOutput(1:ns,2,ipc,ip,el) + cs = cs + ns + + case (rho_dot_ann_ath_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotAthermalAnnihilationOutput(1:ns,1,ipc,ip,el) & + + rhoDotAthermalAnnihilationOutput(1:ns,2,ipc,ip,el) + cs = cs + ns + + case (rho_dot_ann_the_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotThermalAnnihilationOutput(1:ns,1,ipc,ip,el) & + + rhoDotThermalAnnihilationOutput(1:ns,2,ipc,ip,el) + cs = cs + ns + + case (rho_dot_ann_the_edge_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotThermalAnnihilationOutput(1:ns,1,ipc,ip,el) + cs = cs + ns + + case (rho_dot_ann_the_screw_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotThermalAnnihilationOutput(1:ns,2,ipc,ip,el) + cs = cs + ns + + case (rho_dot_edgejogs_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotEdgeJogsOutput(1:ns,ipc,ip,el) + cs = cs + ns + + case (rho_dot_flux_mobile_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotFluxOutput(1:ns,1:4,ipc,ip,el),2) + cs = cs + ns + + case (rho_dot_flux_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotFluxOutput(1:ns,1:4,ipc,ip,el),2) & + + sum(rhoDotFluxOutput(1:ns,5:8,ipc,ip,el)*sign(1.0_pReal,rhoSgl(1:ns,5:8)),2) + cs = cs + ns + + case (rho_dot_flux_edge_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotFluxOutput(1:ns,1:2,ipc,ip,el),2) & + + sum(rhoDotFluxOutput(1:ns,5:6,ipc,ip,el)*sign(1.0_pReal,rhoSgl(1:ns,5:6)),2) + cs = cs + ns + + case (rho_dot_flux_screw_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotFluxOutput(1:ns,3:4,ipc,ip,el),2) & + + sum(rhoDotFluxOutput(1:ns,7:8,ipc,ip,el)*sign(1.0_pReal,rhoSgl(1:ns,7:8)),2) + cs = cs + ns + + case (velocity_edge_pos_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = v(1:ns,1) + cs = cs + ns + + case (velocity_edge_neg_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = v(1:ns,2) + cs = cs + ns + + case (velocity_screw_pos_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = v(1:ns,3) + cs = cs + ns + + case (velocity_screw_neg_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = v(1:ns,4) + cs = cs + ns + + case (slipdirectionx_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = m_currentconf(1,1:ns,1) + cs = cs + ns + + case (slipdirectiony_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = m_currentconf(2,1:ns,1) + cs = cs + ns + + case (slipdirectionz_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = m_currentconf(3,1:ns,1) + cs = cs + ns + + case (slipnormalx_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = n_currentconf(1,1:ns) + cs = cs + ns + + case (slipnormaly_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = n_currentconf(2,1:ns) + cs = cs + ns + + case (slipnormalz_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = n_currentconf(3,1:ns) + cs = cs + ns + + case (fluxdensity_edge_posx_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,1) * v(1:ns,1) * m_currentconf(1,1:ns,1) + cs = cs + ns + + case (fluxdensity_edge_posy_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,1) * v(1:ns,1) * m_currentconf(2,1:ns,1) + cs = cs + ns + + case (fluxdensity_edge_posz_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,1) * v(1:ns,1) * m_currentconf(3,1:ns,1) + cs = cs + ns + + case (fluxdensity_edge_negx_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = - rhoSgl(1:ns,2) * v(1:ns,2) * m_currentconf(1,1:ns,1) + cs = cs + ns + + case (fluxdensity_edge_negy_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = - rhoSgl(1:ns,2) * v(1:ns,2) * m_currentconf(2,1:ns,1) + cs = cs + ns + + case (fluxdensity_edge_negz_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = - rhoSgl(1:ns,2) * v(1:ns,2) * m_currentconf(3,1:ns,1) + cs = cs + ns + + case (fluxdensity_screw_posx_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,3) * v(1:ns,3) * m_currentconf(1,1:ns,2) + cs = cs + ns + + case (fluxdensity_screw_posy_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,3) * v(1:ns,3) * m_currentconf(2,1:ns,2) + cs = cs + ns + + case (fluxdensity_screw_posz_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,3) * v(1:ns,3) * m_currentconf(3,1:ns,2) + cs = cs + ns + + case (fluxdensity_screw_negx_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = - rhoSgl(1:ns,4) * v(1:ns,4) * m_currentconf(1,1:ns,2) + cs = cs + ns + + case (fluxdensity_screw_negy_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = - rhoSgl(1:ns,4) * v(1:ns,4) * m_currentconf(2,1:ns,2) + cs = cs + ns + + case (fluxdensity_screw_negz_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = - rhoSgl(1:ns,4) * v(1:ns,4) * m_currentconf(3,1:ns,2) + cs = cs + ns + + case (maximumdipoleheight_edge_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = dUpper(1:ns,1) + cs = cs + ns + + case (maximumdipoleheight_screw_ID) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = dUpper(1:ns,2) + cs = cs + ns + + case(dislocationstress_ID) + sigma = constitutive_nonlocal_dislocationstress(mapping, Fe, ipc, ip, el) + constitutive_nonlocal_postResults(cs+1_pInt) = sigma(1,1) + constitutive_nonlocal_postResults(cs+2_pInt) = sigma(2,2) + constitutive_nonlocal_postResults(cs+3_pInt) = sigma(3,3) + constitutive_nonlocal_postResults(cs+4_pInt) = sigma(1,2) + constitutive_nonlocal_postResults(cs+5_pInt) = sigma(2,3) + constitutive_nonlocal_postResults(cs+6_pInt) = sigma(3,1) + cs = cs + 6_pInt + + case(accumulatedshear_ID) +! constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = state(ipc,ip,el)%p(iGamma(1:ns,instance)) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = & + plasticState(mapping(2,ipc,ip,el))%state(iGamma(1:ns,instance),mapping(1,ipc,ip,el)) + cs = cs + ns + + end select +enddo outputsLoop + +end function constitutive_nonlocal_postResults + +#else !-------------------------------------------------------------------------------------------------- !> @brief return array of constitutive results !-------------------------------------------------------------------------------------------------- @@ -3781,5 +6208,5 @@ outputsLoop: do o = 1_pInt,phase_Noutput(material_phase(ipc,ip,el)) enddo outputsLoop end function constitutive_nonlocal_postResults - +#endif end module constitutive_nonlocal diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 85df4073a..8c25d70e9 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -1287,7 +1287,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity) case(5_pInt) call crystallite_integrateStateRKCK45() end select - !why not OMP? elementLooping8: do e = FEsolving_execElem(1),FEsolving_execElem(2) myNgrains = homogenization_Ngrains(mesh_element(3,e)) @@ -2409,9 +2408,9 @@ subroutine crystallite_integrateStateAdaptiveEuler() mySizeDotState = plasticState(mappingConstitutive(2,g,i,e))%sizeDotState stateResiduum(1:mySizeDotState,g,i,e) = - 0.5_pReal * plasticState(mappingConstitutive(2,g,i,e))% & dotstate(1:mySizeDotState,mappingConstitutive(1,g,i,e)) * crystallite_subdt(g,i,e) ! contribution to absolute residuum in state - plasticState(mappingConstitutive(2,g,i,e))%state(:,mappingConstitutive(1,g,i,e)) = & - plasticState(mappingConstitutive(2,g,i,e))%state(:,mappingConstitutive(1,g,i,e)) & - + plasticState(mappingConstitutive(2,g,i,e))%dotstate(:,mappingConstitutive(1,g,i,e)) & + plasticState(mappingConstitutive(2,g,i,e))%state(1:mySizeDotState,mappingConstitutive(1,g,i,e)) = & + plasticState(mappingConstitutive(2,g,i,e))%state(1:mySizeDotState,mappingConstitutive(1,g,i,e)) & + + plasticState(mappingConstitutive(2,g,i,e))%dotstate(1:mySizeDotState,mappingConstitutive(1,g,i,e)) & * crystallite_subdt(g,i,e) #endif endif @@ -3335,7 +3334,6 @@ real(pReal), dimension(constitutive_maxSizeDotState) :: & enddo elemLoop enddo crystalliteLooping - end subroutine crystallite_integrateStateFPI