new state related changes

This commit is contained in:
Luv Sharma 2014-06-11 12:27:41 +00:00
parent a54bb6ab24
commit 401b31c951
3 changed files with 421 additions and 152 deletions

View File

@ -137,9 +137,9 @@ subroutine constitutive_init
use constitutive_j2
use constitutive_phenopowerlaw
use constitutive_dislotwin
#ifndef NEWSTATE
use constitutive_titanmod
use constitutive_nonlocal
#ifndef NEWSTATE
use constitutive_nonlocal
#endif
implicit none
integer(pInt), parameter :: FILEUNIT = 200_pInt
@ -210,11 +210,11 @@ subroutine constitutive_init
outputName = PLASTICITY_DISLOTWIN_label
thisOutput => constitutive_dislotwin_output
thisSize => constitutive_dislotwin_sizePostResult
#ifndef NEWSTATE
case (PLASTICITY_TITANMOD_ID)
outputName = PLASTICITY_TITANMOD_label
thisOutput => constitutive_titanmod_output
thisSize => constitutive_titanmod_sizePostResult
#ifndef NEWSTATE
case (PLASTICITY_NONLOCAL_ID)
outputName = PLASTICITY_NONLOCAL_label
thisOutput => constitutive_nonlocal_output
@ -344,25 +344,25 @@ subroutine constitutive_init
case (PLASTICITY_PHENOPOWERLAW_ID)
#ifndef NEWSTATE
allocate(constitutive_state0(g,i,e)%p(constitutive_phenopowerlaw_sizeState(instance)))
allocate(constitutive_partionedState0(g,i,e)%p(constitutive_phenopowerlaw_sizeState(instance)))
allocate(constitutive_subState0(g,i,e)%p(constitutive_phenopowerlaw_sizeState(instance)))
allocate(constitutive_state(g,i,e)%p(constitutive_phenopowerlaw_sizeState(instance)))
allocate(constitutive_state_backup(g,i,e)%p(constitutive_phenopowerlaw_sizeState(instance)))
allocate(constitutive_aTolState(g,i,e)%p(constitutive_phenopowerlaw_sizeState(instance)))
allocate(constitutive_dotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(instance)))
allocate(constitutive_deltaState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(instance)))
allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(instance)))
allocate(constitutive_state0(g,i,e)%p(constitutive_phenopowerlaw_sizeState(instance)),source=0.0_pReal)
allocate(constitutive_partionedState0(g,i,e)%p(constitutive_phenopowerlaw_sizeState(instance)),source=0.0_pReal)
allocate(constitutive_subState0(g,i,e)%p(constitutive_phenopowerlaw_sizeState(instance)),source=0.0_pReal)
allocate(constitutive_state(g,i,e)%p(constitutive_phenopowerlaw_sizeState(instance)),source=0.0_pReal)
allocate(constitutive_state_backup(g,i,e)%p(constitutive_phenopowerlaw_sizeState(instance)),source=0.0_pReal)
allocate(constitutive_aTolState(g,i,e)%p(constitutive_phenopowerlaw_sizeState(instance)),source=0.0_pReal)
allocate(constitutive_dotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(instance)),source=0.0_pReal)
allocate(constitutive_deltaState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(instance)),source=0.0_pReal)
allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(instance)),source=0.0_pReal)
if (any(numerics_integrator == 1_pInt)) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(instance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(instance)))
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(instance)),source=0.0_pReal)
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(instance)),source=0.0_pReal)
endif
if (any(numerics_integrator == 4_pInt)) then
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(instance)))
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(instance)),source=0.0_pReal)
endif
if (any(numerics_integrator == 5_pInt)) then
do s = 1_pInt,6_pInt
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(instance)))
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(instance)),source=0.0_pReal)
enddo
endif
constitutive_state0(g,i,e)%p = constitutive_phenopowerlaw_stateInit(instance)
@ -371,36 +371,38 @@ subroutine constitutive_init
constitutive_sizeDotState(g,i,e) = constitutive_phenopowerlaw_sizeDotState(instance)
#endif
constitutive_sizePostResults(g,i,e) = constitutive_phenopowerlaw_sizePostResults(instance)
#ifndef NEWSTATE
case (PLASTICITY_DISLOTWIN_ID)
allocate(constitutive_state0(g,i,e)%p(constitutive_dislotwin_sizeState(instance)))
allocate(constitutive_partionedState0(g,i,e)%p(constitutive_dislotwin_sizeState(instance)))
allocate(constitutive_subState0(g,i,e)%p(constitutive_dislotwin_sizeState(instance)))
allocate(constitutive_state(g,i,e)%p(constitutive_dislotwin_sizeState(instance)))
allocate(constitutive_state_backup(g,i,e)%p(constitutive_dislotwin_sizeState(instance)))
allocate(constitutive_aTolState(g,i,e)%p(constitutive_dislotwin_sizeState(instance)))
allocate(constitutive_dotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(instance)))
allocate(constitutive_deltaState(g,i,e)%p(constitutive_dislotwin_sizeDotState(instance)))
allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_dislotwin_sizeDotState(instance)))
#ifndef NEWSTATE
allocate(constitutive_state0(g,i,e)%p(constitutive_dislotwin_sizeState(instance)),source=0.0_pReal)
allocate(constitutive_partionedState0(g,i,e)%p(constitutive_dislotwin_sizeState(instance)),source=0.0_pReal)
allocate(constitutive_subState0(g,i,e)%p(constitutive_dislotwin_sizeState(instance)),source=0.0_pReal)
allocate(constitutive_state(g,i,e)%p(constitutive_dislotwin_sizeState(instance)),source=0.0_pReal)
allocate(constitutive_state_backup(g,i,e)%p(constitutive_dislotwin_sizeState(instance)),source=0.0_pReal)
allocate(constitutive_aTolState(g,i,e)%p(constitutive_dislotwin_sizeState(instance)),source=0.0_pReal)
allocate(constitutive_dotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(instance)),source=0.0_pReal)
allocate(constitutive_deltaState(g,i,e)%p(constitutive_dislotwin_sizeDotState(instance)),source=0.0_pReal)
allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_dislotwin_sizeDotState(instance)),source=0.0_pReal)
if (any(numerics_integrator == 1_pInt)) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(instance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_dislotwin_sizeDotState(instance)))
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(instance)),source=0.0_pReal)
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_dislotwin_sizeDotState(instance)),source=0.0_pReal)
endif
if (any(numerics_integrator == 4_pInt)) then
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(instance)))
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(instance)),source=0.0_pReal)
endif
if (any(numerics_integrator == 5_pInt)) then
do s = 1_pInt,6_pInt
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_dislotwin_sizeDotState(instance)))
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_dislotwin_sizeDotState(instance)),source=0.0_pReal)
enddo
endif
constitutive_state0(g,i,e)%p = constitutive_dislotwin_stateInit(instance,material_phase(g,i,e))
constitutive_aTolState(g,i,e)%p = constitutive_dislotwin_aTolState(instance)
constitutive_sizeState(g,i,e) = constitutive_dislotwin_sizeState(instance)
constitutive_sizeDotState(g,i,e) = constitutive_dislotwin_sizeDotState(instance)
#endif
constitutive_sizePostResults(g,i,e) = constitutive_dislotwin_sizePostResults(instance)
case (PLASTICITY_TITANMOD_ID)
#ifndef NEWSTATE
allocate(constitutive_state0(g,i,e)%p(constitutive_titanmod_sizeState(instance)))
allocate(constitutive_partionedState0(g,i,e)%p(constitutive_titanmod_sizeState(instance)))
allocate(constitutive_subState0(g,i,e)%p(constitutive_titanmod_sizeState(instance)))
@ -426,8 +428,10 @@ subroutine constitutive_init
constitutive_aTolState(g,i,e)%p = constitutive_titanmod_aTolState(instance)
constitutive_sizeState(g,i,e) = constitutive_titanmod_sizeState(instance)
constitutive_sizeDotState(g,i,e) = constitutive_titanmod_sizeDotState(instance)
#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)
allocate(constitutive_state0(g,i,e)%p(constitutive_nonlocal_sizeState(instance)))
@ -456,6 +460,7 @@ subroutine constitutive_init
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)
end select
enddo GrainLoop
enddo IPloop
@ -469,7 +474,6 @@ subroutine constitutive_init
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
enddo
#endif
@ -535,10 +539,8 @@ pure function constitutive_homogenizedC(ipc,ip,el)
plasticState,&
#endif
PLASTICITY_DISLOTWIN_ID
#ifndef NEWSTATE
use constitutive_titanmod, only: &
constitutive_titanmod_homogenizedC
#endif
use constitutive_dislotwin, only: &
constitutive_dislotwin_homogenizedC
use lattice, only: &
@ -563,8 +565,13 @@ pure function constitutive_homogenizedC(ipc,ip,el)
constitutive_homogenizedC = constitutive_dislotwin_homogenizedC &
(constitutive_state(ipc,ip,el),ipc,ip,el)
#endif
#ifndef NEWSTATE
case (PLASTICITY_TITANMOD_ID)
#ifdef NEWSTATE
constitutive_homogenizedC = constitutive_titanmod_homogenizedC &
(plasticState(mappingConstitutive(2,ipc,ip,el))%state(:,mappingConstitutive(1,ipc,ip,el)), &
ipc,ip,el)
#else
constitutive_homogenizedC = constitutive_titanmod_homogenizedC(constitutive_state(ipc,ip,el), &
ipc,ip,el)
#endif
@ -589,9 +596,9 @@ subroutine constitutive_microstructure(temperature, Fe, Fp, ipc, ip, el)
#endif
PLASTICITY_TITANMOD_ID, &
PLASTICITY_NONLOCAL_ID
#ifndef NEWSTATE
use constitutive_titanmod, only: &
constitutive_titanmod_microstructure
#ifndef NEWSTATE
use constitutive_nonlocal, only: &
constitutive_nonlocal_microstructure
#endif
@ -621,10 +628,17 @@ subroutine constitutive_microstructure(temperature, Fe, Fp, ipc, ip, el)
call constitutive_dislotwin_microstructure(temperature,constitutive_state(ipc,ip,el), &
ipc,ip,el)
#endif
#ifndef NEWSTATE
case (PLASTICITY_TITANMOD_ID)
#ifdef NEWSTATE
call constitutive_titanmod_microstructure(temperature, &
plasticState(mappingConstitutive(2,ipc,ip,el))%state(:,mappingConstitutive(1,ipc,ip,el)), &
ipc,ip,el)
#else
call constitutive_titanmod_microstructure(temperature,constitutive_state(ipc,ip,el), &
ipc,ip,el)
#endif
#ifndef NEWSTATE
case (PLASTICITY_NONLOCAL_ID)
call constitutive_nonlocal_microstructure(constitutive_state,Fe,Fp,ipc,ip,el)
#endif
@ -657,9 +671,9 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, temperature, ip
constitutive_phenopowerlaw_LpAndItsTangent
use constitutive_dislotwin, only: &
constitutive_dislotwin_LpAndItsTangent
#ifndef NEWSTATE
use constitutive_titanmod, only: &
constitutive_titanmod_LpAndItsTangent
#ifndef NEWSTATE
use constitutive_nonlocal, only: &
constitutive_nonlocal_LpAndItsTangent
#endif
@ -708,10 +722,17 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, temperature, ip
call constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, &
temperature,constitutive_state(ipc,ip,el),ipc,ip,el)
#endif
#ifndef NEWSTATE
case (PLASTICITY_TITANMOD_ID)
#ifdef NEWSTATE
call constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,temperature, &
plasticState(mappingConstitutive(2,ipc,ip,el))%state(:,mappingConstitutive(1,ipc,ip,el)), &
ipc,ip,el)
#else
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)
call constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, &
temperature, constitutive_state(ipc,ip,el), ipc,ip,el)
@ -824,9 +845,9 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, Temperature,
constitutive_phenopowerlaw_dotState
use constitutive_dislotwin, only: &
constitutive_dislotwin_dotState
#ifndef NEWSTATE
use constitutive_titanmod, only: &
constitutive_titanmod_dotState
#ifndef NEWSTATE
use constitutive_nonlocal, only: &
constitutive_nonlocal_dotState
#endif
@ -884,10 +905,17 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, Temperature,
constitutive_dotState(ipc,ip,el)%p = constitutive_dislotwin_dotState(Tstar_v,Temperature,&
constitutive_state(ipc,ip,el), ipc,ip,el)
#endif
#ifndef NEWSTATE
case (PLASTICITY_TITANMOD_ID)
#ifdef NEWSTATE
plasticState(mappingConstitutive(2,ipc,ip,el))%dotState(:,mappingConstitutive(1,ipc,ip,el)) &
= constitutive_titanmod_dotState(Tstar_v,Temperature,&
plasticState(mappingConstitutive(2,ipc,ip,el))% &
state(:,mappingConstitutive(1,ipc,ip,el)), ipc,ip,el)
#else
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)
constitutive_dotState(ipc,ip,el)%p = constitutive_nonlocal_dotState(Tstar_v, FeArray, FpArray, &
Temperature, constitutive_state, constitutive_state0, subdt, &
@ -997,9 +1025,9 @@ function constitutive_postResults(Tstar_v, FeArray, temperature, ipc, ip, el)
constitutive_phenopowerlaw_postResults
use constitutive_dislotwin, only: &
constitutive_dislotwin_postResults
#ifndef NEWSTATE
use constitutive_titanmod, only: &
constitutive_titanmod_postResults
#ifndef NEWSTATE
use constitutive_nonlocal, only: &
constitutive_nonlocal_postResults
#endif
@ -1022,16 +1050,16 @@ function constitutive_postResults(Tstar_v, FeArray, temperature, ipc, ip, el)
select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_NONE_ID)
#ifndef NEWSTATE
case (PLASTICITY_TITANMOD_ID)
#ifdef NEWSTATE
constitutive_postResults = constitutive_titanmod_postResults(&
plasticState(mappingConstitutive(2,ipc,ip,el))% &
state(:,mappingConstitutive(1,ipc,ip,el)),ipc,ip,el)
#else
constitutive_postResults = constitutive_titanmod_postResults(&
constitutive_state(ipc,ip,el),ipc,ip,el)
#endif
case (PLASTICITY_J2_ID)
#ifdef HDF
call constitutive_j2_postResults2(Tstar_v,constitutive_state(ipc,ip,el),ipc,ip,el,1)
#endif
#ifdef NEWSTATE
constitutive_postResults= constitutive_j2_postResults(Tstar_v, &
plasticState(mappingConstitutive(2,ipc,ip,el))% &

View File

@ -219,8 +219,15 @@ subroutine constitutive_titanmod_init(fileUnit)
phase_Noutput, &
PLASTICITY_TITANMOD_label, &
PLASTICITY_TITANMOD_ID, &
#ifdef NEWSTATE
plasticState, &
#endif
MATERIAL_partPhase
use lattice
#ifdef NEWSTATE
use numerics,only: &
numerics_integrator
#endif
implicit none
integer(pInt), intent(in) :: fileUnit
@ -238,6 +245,10 @@ subroutine constitutive_titanmod_init(fileUnit)
Nchunks_SlipFamilies, Nchunks_TwinFamilies, &
mySize, &
maxTotalNslip,maxTotalNtwin, maxNinstance
#ifdef NEWSTATE
integer(pInt) :: sizeState, sizeDotState
#endif
integer(pInt) :: NofMyPhase
character(len=65536) :: &
tag = '', &
line = ''
@ -734,6 +745,16 @@ subroutine constitutive_titanmod_init(fileUnit)
constitutive_titanmod_sizeDotState(instance)+ &
size(constitutive_titanmod_listDependentSlipStates)*ns + &
size(constitutive_titanmod_listDependentTwinStates)*nt
#ifdef NEWSTATE
sizeDotState = &
size(constitutive_titanmod_listBasicSlipStates)*ns + &
size(constitutive_titanmod_listBasicTwinStates)*nt
sizeState = &
constitutive_titanmod_sizeDotState(instance)+ &
size(constitutive_titanmod_listDependentSlipStates)*ns + &
size(constitutive_titanmod_listDependentTwinStates)*nt
#endif
!--------------------------------------------------------------------------------------------------
! determine size of postResults array
@ -767,7 +788,29 @@ subroutine constitutive_titanmod_init(fileUnit)
constitutive_titanmod_sizePostResults(instance) = constitutive_titanmod_sizePostResults(instance) + mySize
endif outputFound
enddo outputsLoop
#ifdef NEWSTATE
! Determine size of state array
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
!--------------------------------------------------------------------------------------------------
! construction of the twin elasticity matrices
do j=1_pInt,lattice_maxNtwinFamily
@ -959,6 +1002,122 @@ end subroutine constitutive_titanmod_init
!--------------------------------------------------------------------------------------------------
!> @brief sets the initial microstructural state for a given instance of this plasticity
!--------------------------------------------------------------------------------------------------
#ifdef NEWSTATE
subroutine constitutive_titanmod_stateInit(instance,phase)
use lattice, only: &
lattice_maxNslipFamily, &
lattice_maxNtwinFamily, &
lattice_mu
use material, only: &
plasticState
implicit none
integer(pInt), intent(in) :: instance !< number specifying the instance of the plasticity
integer(pInt), intent(in) :: phase !< number specifying the phase of the plasticity
integer(pInt) :: &
s,s0,s1, &
t,t0,t1, &
ns,nt,f
real(pReal), dimension(constitutive_titanmod_totalNslip(instance)) :: &
rho_edge0, &
rho_screw0, &
shear_system0, &
segment_edge0, &
segment_screw0, &
resistance_edge0, &
resistance_screw0
real(pReal), dimension(constitutive_titanmod_totalNtwin(instance)) :: &
twingamma_dot0, &
resistance_twin0
real(pReal), dimension(plasticState(phase)%sizeState) :: tempState
ns = constitutive_titanmod_totalNslip(instance)
nt = constitutive_titanmod_totalNtwin(instance)
tempState = 0.0_pReal
!--------------------------------------------------------------------------------------------------
! initialize basic slip state variables for slip
s1 = 0_pInt
do f = 1_pInt,lattice_maxNslipFamily
s0 = s1 + 1_pInt
s1 = s0 + constitutive_titanmod_Nslip(f,instance) - 1_pInt
do s = s0,s1
rho_edge0(s) = constitutive_titanmod_rho_edge0(f,instance)
rho_screw0(s) = constitutive_titanmod_rho_screw0(f,instance)
shear_system0(s) = 0.0_pReal
enddo
enddo
!--------------------------------------------------------------------------------------------------
! initialize basic slip state variables for twin
t1 = 0_pInt
do f = 1_pInt,lattice_maxNtwinFamily
t0 = t1 + 1_pInt
t1 = t0 + constitutive_titanmod_Ntwin(f,instance) - 1_pInt
do t = t0,t1
twingamma_dot0(t)=0.0_pReal
enddo
enddo
!--------------------------------------------------------------------------------------------------
! initialize dependent slip microstructural variables
forall (s = 1_pInt:ns)
segment_edge0(s) = constitutive_titanmod_CeLambdaSlipPerSlipSys(s,instance)/ &
sqrt(dot_product((rho_edge0),constitutive_titanmod_forestProjectionEdge(1:ns,s,instance))+ &
dot_product((rho_screw0),constitutive_titanmod_forestProjectionScrew(1:ns,s,instance)))
segment_screw0(s) = constitutive_titanmod_CsLambdaSlipPerSlipSys(s,instance)/ &
sqrt(dot_product((rho_edge0),constitutive_titanmod_forestProjectionEdge(1:ns,s,instance))+ &
dot_product((rho_screw0),constitutive_titanmod_forestProjectionScrew(1:ns,s,instance)))
resistance_edge0(s) = &
lattice_mu(phase)*constitutive_titanmod_burgersPerSlipSys(s,instance)* &
sqrt(dot_product((rho_edge0),constitutive_titanmod_interactionMatrix_ee(1:ns,s,instance))+ &
dot_product((rho_screw0),constitutive_titanmod_interactionMatrix_es(1:ns,s,instance)))
resistance_screw0(s) = &
lattice_mu(phase)*constitutive_titanmod_burgersPerSlipSys(s,instance)* &
sqrt(dot_product((rho_edge0),constitutive_titanmod_interactionMatrix_es(1:ns,s,instance))+ &
dot_product((rho_screw0), constitutive_titanmod_interactionMatrix_ss(1:ns,s,instance)))
end forall
forall (t = 1_pInt:nt) &
resistance_twin0(t) = 0.0_pReal
tempState = 0.0_pReal
tempState (1:ns) = rho_edge0
tempState (1_pInt*ns+1_pInt:2_pInt*ns) = rho_screw0
tempState (2_pInt*ns+1_pInt:3_pInt*ns) = shear_system0
tempState (3_pInt*ns+1_pInt:3_pInt*ns+nt) = twingamma_dot0
tempState (3_pInt*ns+nt+1_pInt:4_pInt*ns+nt) = segment_edge0
tempState (4_pInt*ns+nt+1_pInt:5_pInt*ns+nt) = segment_screw0
tempState (5_pInt*ns+nt+1_pInt:6_pInt*ns+nt) = resistance_edge0
tempState (6_pInt*ns+nt+1_pInt:7_pInt*ns+nt) = resistance_screw0
tempState (7_pInt*ns+nt+1_pInt:7_pInt*ns+2_pInt*nt)=resistance_twin0
end subroutine constitutive_titanmod_stateInit
!--------------------------------------------------------------------------------------------------
!> @brief sets the relevant state values for a given instance of this plasticity
!--------------------------------------------------------------------------------------------------
subroutine constitutive_titanmod_aTolState(phase,instance)
use material, only: &
plasticState
implicit none
integer(pInt), intent(in) :: &
instance, &
phase
! real(pReal), dimension(size(plasticState(phase)%aTolState(:))) :: tempTol
real(pReal), dimension(plasticState(phase)%sizeState) :: tempTol
tempTol = 0.0_pReal
tempTol = constitutive_titanmod_aTolRho(instance)
end subroutine constitutive_titanmod_aTolState
#else
pure function constitutive_titanmod_stateInit(instance,phase)
use lattice, only: &
lattice_maxNslipFamily, &
@ -1066,7 +1225,7 @@ pure function constitutive_titanmod_aTolState(instance)
end function constitutive_titanmod_aTolState
#endif
!--------------------------------------------------------------------------------------------------
!> @brief returns the homogenized elasticity matrix
!--------------------------------------------------------------------------------------------------
@ -1090,8 +1249,17 @@ implicit none
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
type(p_vec), intent(in) :: &
#ifdef NEWSTATE
real(pReal), dimension(:), intent(in) :: &
state
real(pReal), dimension(size(state)) :: &
tempState
#else
type(p_vec), intent(in) :: &
state !< microstructure state
real(pReal), dimension(size(state%p)) :: &
tempState
#endif
real(pReal), dimension(constitutive_titanmod_totalNtwin(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
volumefraction_PerTwinSys
integer(pInt) :: &
@ -1101,7 +1269,13 @@ real(pReal), dimension(constitutive_titanmod_totalNtwin(phase_plasticityInstance
i
real(pReal) :: &
sumf
tempState = 0.0_pReal
#ifdef NEWSTATE
tempState=state
#else
tempState = state%p
#endif
!--------------------------------------------------------------------------------------------------
! shortened notation
phase = material_phase(ipc,ip,el)
@ -1112,7 +1286,7 @@ real(pReal), dimension(constitutive_titanmod_totalNtwin(phase_plasticityInstance
!--------------------------------------------------------------------------------------------------
! total twin volume fraction
do i=1_pInt,nt
volumefraction_PerTwinSys(i)=state%p(3_pInt*ns+i)/ &
volumefraction_PerTwinSys(i)=tempState(3_pInt*ns+i)/ &
constitutive_titanmod_twinshearconstant_PerTwinSys(i,instance)
enddo
sumf = sum(abs(volumefraction_PerTwinSys(1:nt))) ! safe for nt == 0
@ -1152,8 +1326,17 @@ subroutine constitutive_titanmod_microstructure(temperature,state,ipc,ip,el)
el !< element
real(pReal), intent(in) :: &
temperature !< temperature at IP
type(p_vec), intent(inout) :: &
#ifdef NEWSTATE
real(pReal), dimension(:), intent(inout) :: &
state
real(pReal), dimension(size(state)) :: &
tempState
#else
type(p_vec), intent(inout) :: &
state !< microstructure state
real(pReal), dimension(size(state%p)) :: &
tempState
#endif
integer(pInt) :: &
instance, &
@ -1166,6 +1349,12 @@ subroutine constitutive_titanmod_microstructure(temperature,state,ipc,ip,el)
volumefraction_PerTwinSys
!--------------------------------------------------------------------------------------------------
#ifdef NEWSTATE
tempState=state
#else
tempState = state%p
#endif
!Shortened notation
phase = material_phase(ipc,ip,el)
instance = phase_plasticityInstance(phase)
@ -1175,7 +1364,7 @@ subroutine constitutive_titanmod_microstructure(temperature,state,ipc,ip,el)
!--------------------------------------------------------------------------------------------------
! total twin volume fraction
forall (i = 1_pInt:nt) &
volumefraction_PerTwinSys(i)=state%p(3_pInt*ns+i)/ &
volumefraction_PerTwinSys(i)=tempState(3_pInt*ns+i)/ &
constitutive_titanmod_twinshearconstant_PerTwinSys(i,instance)
sumf = sum(abs(volumefraction_PerTwinSys(1:nt))) ! safe for nt == 0
@ -1185,44 +1374,50 @@ subroutine constitutive_titanmod_microstructure(temperature,state,ipc,ip,el)
!--------------------------------------------------------------------------------------------------
! average segment length for edge dislocations in matrix
forall (s = 1_pInt:ns) &
state%p(3_pInt*ns+nt+s) = constitutive_titanmod_CeLambdaSlipPerSlipSys(s,instance)/ &
sqrt(dot_product(state%p(1:ns), &
tempState(3_pInt*ns+nt+s) = constitutive_titanmod_CeLambdaSlipPerSlipSys(s,instance)/ &
sqrt(dot_product(tempState(1:ns), &
constitutive_titanmod_forestProjectionEdge(1:ns,s,instance))+ &
dot_product(state%p(ns+1_pInt:2_pInt*ns), &
dot_product(tempState(ns+1_pInt:2_pInt*ns), &
constitutive_titanmod_forestProjectionScrew(1:ns,s,instance)))
!--------------------------------------------------------------------------------------------------
! average segment length for screw dislocations in matrix
forall (s = 1_pInt:ns) &
state%p(4_pInt*ns+nt+s) = constitutive_titanmod_CsLambdaSlipPerSlipSys(s,instance)/ &
sqrt(dot_product(state%p(1:ns), &
tempState(4_pInt*ns+nt+s) = constitutive_titanmod_CsLambdaSlipPerSlipSys(s,instance)/ &
sqrt(dot_product(tempState(1:ns), &
constitutive_titanmod_forestProjectionEdge(1:ns,s,instance))+ &
dot_product(state%p(ns+1_pInt:2_pInt*ns), &
dot_product(tempState(ns+1_pInt:2_pInt*ns), &
constitutive_titanmod_forestProjectionScrew(1:ns,s,instance)))
!--------------------------------------------------------------------------------------------------
! threshold stress or slip resistance for edge dislocation motion
forall (s = 1_pInt:ns) &
state%p(5_pInt*ns+nt+s) = &
tempState(5_pInt*ns+nt+s) = &
lattice_mu(phase)*constitutive_titanmod_burgersPerSlipSys(s,instance)*&
sqrt(dot_product((state%p(1:ns)),&
sqrt(dot_product((tempState(1:ns)),&
constitutive_titanmod_interactionMatrix_ee(1:ns,s,instance))+ &
dot_product((state%p(ns+1_pInt:2_pInt*ns)),&
dot_product((tempState(ns+1_pInt:2_pInt*ns)),&
constitutive_titanmod_interactionMatrix_es(1:ns,s,instance)))
!--------------------------------------------------------------------------------------------------
! threshold stress or slip resistance for screw dislocation motion
forall (s = 1_pInt:ns) &
state%p(6_pInt*ns+nt+s) = &
tempState(6_pInt*ns+nt+s) = &
lattice_mu(phase)*constitutive_titanmod_burgersPerSlipSys(s,instance)*&
sqrt(dot_product((state%p(1:ns)),&
sqrt(dot_product((tempState(1:ns)),&
constitutive_titanmod_interactionMatrix_es(1:ns,s,instance))+ &
dot_product((state%p(ns+1_pInt:2_pInt*ns)),&
dot_product((tempState(ns+1_pInt:2_pInt*ns)),&
constitutive_titanmod_interactionMatrix_ss(1:ns,s,instance)))
!--------------------------------------------------------------------------------------------------
! threshold stress or slip resistance for dislocation motion in twin
forall (t = 1_pInt:nt) &
state%p(7_pInt*ns+nt+t) = &
tempState(7_pInt*ns+nt+t) = &
lattice_mu(phase)*constitutive_titanmod_burgersPerTwinSys(t,instance)*&
(dot_product((abs(state%p(2_pInt*ns+1_pInt:2_pInt*ns+nt))),&
(dot_product((abs(tempState(2_pInt*ns+1_pInt:2_pInt*ns+nt))),&
constitutive_titanmod_interactionMatrixTwinTwin(1:nt,t,instance)))
#ifdef NEWSTATE
state=tempState
#else
state%p = tempState
#endif
end subroutine constitutive_titanmod_microstructure
@ -1270,8 +1465,17 @@ subroutine constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,&
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
type(p_vec), intent(inout) :: &
#ifdef NEWSTATE
real(pReal), dimension(:), intent(inout) :: &
state
real(pReal), dimension(size(state)) :: &
tempState
#else
type(p_vec), intent(inout) :: &
state !< microstructure state
real(pReal), dimension(size(state%p)) :: &
tempState
#endif
integer(pInt) :: &
index_myFamily, instance,phase, &
ns,nt, &
@ -1288,6 +1492,14 @@ subroutine constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,&
gdot_slip_edge, gdot_slip_screw
real(pReal), dimension(constitutive_titanmod_totalNtwin(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
gdot_twin,dgdot_dtautwin,tau_twin, volumefraction_PerTwinSys
#ifdef NEWSTATE
tempState=state
#else
tempState = state%p
#endif
!--------------------------------------------------------------------------------------------------
! shortened notation
@ -1297,7 +1509,7 @@ subroutine constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,&
nt = constitutive_titanmod_totalNtwin(instance)
do i=1_pInt,nt
volumefraction_PerTwinSys(i)=state%p(3_pInt*ns+i)/ &
volumefraction_PerTwinSys(i)=tempState(3_pInt*ns+i)/ &
constitutive_titanmod_twinshearconstant_PerTwinSys(i,instance)
enddo
@ -1325,7 +1537,7 @@ subroutine constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,&
tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,phase))
if(lattice_structure(phase)==LATTICE_hex_ID) then ! only for prismatic and pyr <a> systems in hex
screwvelocity_prefactor=constitutive_titanmod_debyefrequency(instance)* &
state%p(4_pInt*ns+nt+j)*(constitutive_titanmod_burgersPerSlipSys(j,instance)/ &
tempState(4_pInt*ns+nt+j)*(constitutive_titanmod_burgersPerSlipSys(j,instance)/ &
constitutive_titanmod_kinkcriticallength_PerSlipSys(j,instance))**2
!* Stress ratio for screw ! No slip resistance for screw dislocations, only Peierls stress
@ -1350,7 +1562,7 @@ subroutine constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,&
else ! if the structure is not hex or the slip family is basal
screwvelocity_prefactor=constitutive_titanmod_v0s_PerSlipSys(j,instance)
bottomstress_screw=constitutive_titanmod_tau0s_PerSlipSys(j,instance)+state%p(6*ns+nt+j)
bottomstress_screw=constitutive_titanmod_tau0s_PerSlipSys(j,instance)+tempState(6*ns+nt+j)
StressRatio_screw_p = ((abs(tau_slip(j)))/( bottomstress_screw ))**constitutive_titanmod_ps_PerSlipSys(j,instance)
if((1.0_pReal-StressRatio_screw_p)>0.001_pReal) then
@ -1368,7 +1580,7 @@ subroutine constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,&
endif
!* Stress ratio for edge
bottomstress_edge=constitutive_titanmod_tau0e_PerSlipSys(j,instance)+state%p(5*ns+nt+j)
bottomstress_edge=constitutive_titanmod_tau0e_PerSlipSys(j,instance)+tempState(5*ns+nt+j)
StressRatio_edge_p = ((abs(tau_slip(j)))/ &
( bottomstress_edge) &
)**constitutive_titanmod_pe_PerSlipSys(j,instance)
@ -1394,29 +1606,29 @@ subroutine constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,&
constitutive_titanmod_qe_PerSlipSys(j,instance))
!* Shear rates due to edge slip
gdot_slip_edge(j) = constitutive_titanmod_burgersPerSlipSys(j,instance)*(state%p(j)* &
gdot_slip_edge(j) = constitutive_titanmod_burgersPerSlipSys(j,instance)*(tempState(j)* &
edge_velocity(j))* sign(1.0_pReal,tau_slip(j))
!* Shear rates due to screw slip
gdot_slip_screw(j) = constitutive_titanmod_burgersPerSlipSys(j,instance)*(state%p(ns+j) * &
gdot_slip_screw(j) = constitutive_titanmod_burgersPerSlipSys(j,instance)*(tempState(ns+j) * &
screw_velocity(j))* sign(1.0_pReal,tau_slip(j))
!Total shear rate
gdot_slip(j) = gdot_slip_edge(j) + gdot_slip_screw(j)
state%p(7*ns+2*nt+j)=edge_velocity(j)
state%p(8*ns+2*nt+j)=screw_velocity(j)
state%p(9*ns+2*nt+j)=tau_slip(j)
state%p(10*ns+2*nt+j)=gdot_slip_edge(j)
state%p(11*ns+2*nt+j)=gdot_slip_screw(j)
state%p(12*ns+2*nt+j)=StressRatio_edge_p
state%p(13*ns+2*nt+j)=StressRatio_screw_p
tempState(7*ns+2*nt+j)=edge_velocity(j)
tempState(8*ns+2*nt+j)=screw_velocity(j)
tempState(9*ns+2*nt+j)=tau_slip(j)
tempState(10*ns+2*nt+j)=gdot_slip_edge(j)
tempState(11*ns+2*nt+j)=gdot_slip_screw(j)
tempState(12*ns+2*nt+j)=StressRatio_edge_p
tempState(13*ns+2*nt+j)=StressRatio_screw_p
!* Derivatives of shear rates
dgdot_dtauslip(j) = constitutive_titanmod_burgersPerSlipSys(j,instance)*(( &
( &
( &
( &
(edge_velocity(j)*state%p(j))) * &
(edge_velocity(j)*tempState(j))) * &
BoltzmannRatioedge*&
constitutive_titanmod_pe_PerSlipSys(j,instance)* &
constitutive_titanmod_qe_PerSlipSys(j,instance) &
@ -1429,7 +1641,7 @@ subroutine constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,&
( &
( &
( &
(state%p(ns+j) * screw_velocity(j)) * &
(tempState(ns+j) * screw_velocity(j)) * &
BoltzmannRatioscrew* &
constitutive_titanmod_ps_PerSlipSys(j,instance)* &
constitutive_titanmod_qs_PerSlipSys(j,instance) &
@ -1484,7 +1696,7 @@ subroutine constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,&
!* Stress ratio for edge
twinStressRatio_p = ((abs(tau_twin(j)))/ &
( constitutive_titanmod_twintau0_PerTwinSys(j,instance)+state%p(7*ns+nt+j)) &
( constitutive_titanmod_twintau0_PerTwinSys(j,instance)+tempState(7*ns+nt+j)) &
)**constitutive_titanmod_twinp_PerTwinSys(j,instance)
if((1.0_pReal-twinStressRatio_p)>0.001_pReal) then
@ -1494,7 +1706,7 @@ subroutine constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,&
endif
twinStressRatio_pminus1 = ((abs(tau_twin(j)))/ &
( constitutive_titanmod_twintau0_PerTwinSys(j,instance)+state%p(7*ns+nt+j)) &
( constitutive_titanmod_twintau0_PerTwinSys(j,instance)+tempState(7*ns+nt+j)) &
)**(constitutive_titanmod_twinp_PerTwinSys(j,instance)-1.0_pReal)
!* Boltzmann ratio
@ -1538,6 +1750,11 @@ subroutine constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,&
enddo twinFamiliesLoop
dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333)
#ifdef NEWSTATE
tempState=state
#else
tempState = state%p
#endif
end subroutine constitutive_titanmod_LpAndItsTangent
@ -1571,10 +1788,21 @@ implicit none
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
#ifdef NEWSTATE
real(pReal), dimension(:), intent(in) :: &
state
real(pReal), dimension(size(state)) :: &
tempState
real(pReal), dimension(size(state)) :: &
constitutive_titanmod_dotState
#else
type(p_vec), intent(in) :: &
state !< microstructure state
real(pReal), dimension(size(state%p)) :: &
tempState
real(pReal), dimension(constitutive_titanmod_sizeDotState(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
constitutive_titanmod_dotState
#endif
integer(pInt) :: &
index_myFamily, instance,phase, &
@ -1601,7 +1829,7 @@ implicit none
nt = constitutive_titanmod_totalNtwin(instance)
do i=1_pInt,nt
volumefraction_PerTwinSys(i)=state%p(3_pInt*ns+i)/ &
volumefraction_PerTwinSys(i)=tempState(3_pInt*ns+i)/ &
constitutive_titanmod_twinshearconstant_PerTwinSys(i,instance)
enddo
@ -1617,13 +1845,13 @@ implicit none
j = j+1_pInt
DotRhoEdgeGeneration(j) = & ! multiplication of edge dislocations
state%p(ns+j)*state%p(8*ns+2*nt+j)/state%p(4*ns+nt+j)
tempState(ns+j)*tempState(8*ns+2*nt+j)/tempState(4*ns+nt+j)
DotRhoScrewGeneration(j) = & ! multiplication of screw dislocations
state%p(j)*state%p(7*ns+2*nt+j)/state%p(3*ns+nt+j)
DotRhoEdgeAnnihilation(j) = -((state%p(j))**2)* & ! annihilation of edge dislocations
constitutive_titanmod_capre_PerSlipSys(j,instance)*state%p(7*ns+2*nt+j)*0.5_pReal
DotRhoScrewAnnihilation(j) = -((state%p(ns+j))**2)* & ! annihilation of screw dislocations
constitutive_titanmod_caprs_PerSlipSys(j,instance)*state%p(8*ns+2*nt+j)*0.5_pReal
tempState(j)*tempState(7*ns+2*nt+j)/tempState(3*ns+nt+j)
DotRhoEdgeAnnihilation(j) = -((tempState(j))**2)* & ! annihilation of edge dislocations
constitutive_titanmod_capre_PerSlipSys(j,instance)*tempState(7*ns+2*nt+j)*0.5_pReal
DotRhoScrewAnnihilation(j) = -((tempState(ns+j))**2)* & ! annihilation of screw dislocations
constitutive_titanmod_caprs_PerSlipSys(j,instance)*tempState(8*ns+2*nt+j)*0.5_pReal
constitutive_titanmod_dotState(j) = & ! edge dislocation density rate of change
DotRhoEdgeGeneration(j)+DotRhoEdgeAnnihilation(j)
@ -1631,7 +1859,7 @@ implicit none
DotRhoScrewGeneration(j)+DotRhoScrewAnnihilation(j)
constitutive_titanmod_dotState(2*ns+j) = & ! sum of shear due to edge and screw
state%p(10*ns+2*nt+j)+state%p(11*ns+2*nt+j)
tempState(10*ns+2*nt+j)+tempState(11*ns+2*nt+j)
enddo
enddo slipFamiliesLoop
@ -1647,7 +1875,7 @@ implicit none
!* Stress ratio for edge
twinStressRatio_p = ((abs(tau_twin(j)))/ &
( constitutive_titanmod_twintau0_PerTwinSys(j,instance)+state%p(7*ns+nt+j)) &
( constitutive_titanmod_twintau0_PerTwinSys(j,instance)+tempState(7*ns+nt+j)) &
)**(constitutive_titanmod_twinp_PerTwinSys(j,instance))
@ -1691,8 +1919,17 @@ pure function constitutive_titanmod_postResults(state,ipc,ip,el)
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
type(p_vec), intent(in) :: &
#ifdef NEWSTATE
real(pReal), dimension(:), intent(in) :: &
state
real(pReal), dimension(size(state)) :: &
tempState
#else
type(p_vec), intent(in) :: &
state !< microstructure state
real(pReal), dimension(size(state%p)) :: &
tempState
#endif
real(pReal), dimension(constitutive_titanmod_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
constitutive_titanmod_postResults
@ -1706,6 +1943,11 @@ pure function constitutive_titanmod_postResults(state,ipc,ip,el)
volumefraction_PerTwinSys
!--------------------------------------------------------------------------------------------------
#ifdef NEWSTATE
tempState=state
#else
tempState = state%p
#endif
! shortened notation
phase = material_phase(ipc,ip,el)
instance = phase_plasticityInstance(phase)
@ -1713,7 +1955,7 @@ pure function constitutive_titanmod_postResults(state,ipc,ip,el)
nt = constitutive_titanmod_totalNtwin(instance)
do i=1_pInt,nt
volumefraction_PerTwinSys(i)=state%p(3_pInt*ns+i)/ &
volumefraction_PerTwinSys(i)=tempState(3_pInt*ns+i)/ &
constitutive_titanmod_twinshearconstant_PerTwinSys(i,instance)
enddo
@ -1728,91 +1970,91 @@ pure function constitutive_titanmod_postResults(state,ipc,ip,el)
do o = 1_pInt,phase_Noutput(material_phase(ipc,ip,el))
select case(constitutive_titanmod_outputID(o,instance))
case (rhoedge_ID)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = state%p(1_pInt:ns)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = tempState(1_pInt:ns)
c = c + ns
case (rhoscrew_ID)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = state%p(ns+1_pInt:2_pInt*ns)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = tempState(ns+1_pInt:2_pInt*ns)
c = c + ns
case (segment_edge_ID)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = state%p(3_pInt*ns+nt+1_pInt:4_pInt*ns+nt)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = tempState(3_pInt*ns+nt+1_pInt:4_pInt*ns+nt)
c = c + ns
case (segment_screw_ID)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = state%p(4_pInt*ns+nt+1_pInt:5_pInt*ns+nt)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = tempState(4_pInt*ns+nt+1_pInt:5_pInt*ns+nt)
c = c + ns
case (resistance_edge_ID)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = state%p(5_pInt*ns+nt+1_pInt:6_pInt*ns+nt)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = tempState(5_pInt*ns+nt+1_pInt:6_pInt*ns+nt)
c = c + ns
case (resistance_screw_ID)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = state%p(6_pInt*ns+nt+1_pInt:7_pInt*ns+nt)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = tempState(6_pInt*ns+nt+1_pInt:7_pInt*ns+nt)
c = c + ns
case (velocity_edge_ID)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = state%p(7*ns+2*nt+1:8*ns+2*nt)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = tempState(7*ns+2*nt+1:8*ns+2*nt)
c = c + ns
case (velocity_screw_ID)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = state%p(8*ns+2*nt+1:9*ns+2*nt)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = tempState(8*ns+2*nt+1:9*ns+2*nt)
c = c + ns
case (tau_slip_ID)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state%p(9*ns+2*nt+1:10*ns+2*nt))
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(tempState(9*ns+2*nt+1:10*ns+2*nt))
c = c + ns
case (gdot_slip_edge_ID)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state%p(10*ns+2*nt+1:11*ns+2*nt))
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(tempState(10*ns+2*nt+1:11*ns+2*nt))
c = c + ns
case (gdot_slip_screw_ID)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state%p(11*ns+2*nt+1:12*ns+2*nt))
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(tempState(11*ns+2*nt+1:12*ns+2*nt))
c = c + ns
case (gdot_slip_ID)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state%p(10*ns+2*nt+1:11*ns+2*nt)) + &
abs(state%p(11*ns+2*nt+1:12*ns+2*nt))
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(tempState(10*ns+2*nt+1:11*ns+2*nt)) + &
abs(tempState(11*ns+2*nt+1:12*ns+2*nt))
c = c + ns
case (stressratio_edge_p_ID)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state%p(12*ns+2*nt+1:13*ns+2*nt))
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(tempState(12*ns+2*nt+1:13*ns+2*nt))
c = c + ns
case (stressratio_screw_p_ID)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state%p(13*ns+2*nt+1:14*ns+2*nt))
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(tempState(13*ns+2*nt+1:14*ns+2*nt))
c = c + ns
case (shear_system_ID)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state%p(2*ns+1:3*ns))
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(tempState(2*ns+1:3*ns))
c = c + ns
case (shear_basal_ID)
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state%p(2*ns+1:2*ns+3)))
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(tempState(2*ns+1:2*ns+3)))
c = c + 1_pInt
case (shear_prism_ID)
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state%p(2*ns+4:2*ns+6)))
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(tempState(2*ns+4:2*ns+6)))
c = c + 1_pInt
case (shear_pyra_ID)
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state%p(2*ns+7:2*ns+12)))
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(tempState(2*ns+7:2*ns+12)))
c = c + 1_pInt
case (shear_pyrca_ID)
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state%p(2*ns+13:2*ns+24)))
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(tempState(2*ns+13:2*ns+24)))
c = c + 1_pInt
case (rhoedge_basal_ID)
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state%p(1:3))
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(tempState(1:3))
c = c + 1_pInt
case (rhoedge_prism_ID)
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state%p(4:6))
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(tempState(4:6))
c = c + 1_pInt
case (rhoedge_pyra_ID)
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state%p(7:12))
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(tempState(7:12))
c = c + 1_pInt
case (rhoedge_pyrca_ID)
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state%p(13:24))
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(tempState(13:24))
c = c + 1_pInt
case (rhoscrew_basal_ID)
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state%p(ns+1:ns+3))
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(tempState(ns+1:ns+3))
c = c + 1_pInt
case (rhoscrew_prism_ID)
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state%p(ns+4:ns+6))
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(tempState(ns+4:ns+6))
c = c + 1_pInt
case (rhoscrew_pyra_ID)
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state%p(ns+7:ns+12))
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(tempState(ns+7:ns+12))
c = c + 1_pInt
case (rhoscrew_pyrca_ID)
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state%p(ns+13:ns+24))
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(tempState(ns+13:ns+24))
c = c + 1_pInt
case (shear_total_ID)
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state%p(2*ns+1:3*ns)))
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(tempState(2*ns+1:3*ns)))
c = c + 1_pInt
case (twin_fraction_ID)
constitutive_titanmod_postResults(c+1_pInt:c+nt) = abs(volumefraction_PerTwinSys(1:nt))

View File

@ -425,7 +425,6 @@ subroutine crystallite_init(temperature)
enddo
enddo
!$OMP END PARALLEL DO
call crystallite_stressAndItsTangent(.true.,.false.) ! request elastic answers
crystallite_fallbackdPdF = crystallite_dPdF ! use initial elastic stiffness as fallback
@ -1749,9 +1748,9 @@ subroutine crystallite_integrateStateRKCK45()
homogenization_maxNgrains
use constitutive, only: &
constitutive_collectDotState, &
constitutive_maxSizeDotState, &
#ifndef NEWSTATE
constitutive_sizeDotState, &
constitutive_maxSizeDotState, &
constitutive_state, &
constitutive_aTolState, &
constitutive_subState0, &
@ -2245,7 +2244,7 @@ subroutine crystallite_integrateStateRKCK45()
!$OMP DO
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then
crystallite_converged(g,i,e) = .true. ! if still "to do" then converged per definitionem
crystallite_converged(g,i,e) = .true. ! if still "to do" then converged per definitionen
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (distributionState)
debug_StateLoopDistribution(6,numerics_integrationMode) = &
@ -2262,7 +2261,7 @@ subroutine crystallite_integrateStateRKCK45()
! --- nonlocal convergence check ---
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) &
write(6,'(a,i8,a,i2,/)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), ' grains converged' ! if not requesting Integration of just a single IP
write(6,'(a,i8,a,i2,/)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), ' grains converged' ! if not requesting Integration of just a single IP
if ((.not. singleRun) .and. any(.not. crystallite_converged .and. .not. crystallite_localPlasticity)) & ! any non-local not yet converged (or broken)...
crystallite_converged = crystallite_converged .and. crystallite_localPlasticity ! ...restart all non-local as not converged
@ -2717,7 +2716,6 @@ eIter = FEsolving_execElem(1:2)
enddo
singleRun = (eIter(1) == eIter(2) .and. iIter(1,eIter(1)) == iIter(2,eIter(2)))
if (numerics_integrationMode == 1_pInt) then
!$OMP PARALLEL
@ -3010,7 +3008,7 @@ real(pReal), dimension(constitutive_maxSizeDotState) :: &
endif
#else
if ( any(plasticState(mappingConstitutive(2,g,i,e))%dotState(:,mappingConstitutive(1,g,i,e))/= &
plasticState(mappingConstitutive(2,g,i,e))%dotState(:,mappingConstitutive(1,g,i,e)))) then ! NaN occured in dotState
plasticState(mappingConstitutive(2,g,i,e))%dotState(:,mappingConstitutive(1,g,i,e)))) then ! NaN occured in dotState
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken is a non-local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals done (and broken)
@ -3036,17 +3034,12 @@ real(pReal), dimension(constitutive_maxSizeDotState) :: &
+ constitutive_dotState(g,i,e)%p(1:mySizeDotState) &
* crystallite_subdt(g,i,e)
#else
! write(6,*) size(plasticState); flush(6)
! write(6,*) g,i,e;flush(6)
! write(6,*) mappingConstitutive(2,g,i,e),mappingConstitutive(1,g,i,e); flush(6)
! write(6,*) size(plasticState(mappingConstitutive(2,g,i,e))%state(:,mappingConstitutive(1,g,i,e)));flush(6)
! write(6,*) size(plasticState(mappingConstitutive(2,g,i,e))%substate0(:,mappingConstitutive(1,g,i,e)));flush(6)
! write(6,*) size(plasticState(mappingConstitutive(2,g,i,e))%dotstate(:,mappingConstitutive(1,g,i,e)));flush(6)
plasticState(mappingConstitutive(2,g,i,e))%state(:,mappingConstitutive(1,g,i,e)) = &
plasticState(mappingConstitutive(2,g,i,e))%subState0(:,mappingConstitutive(1,g,i,e)) &
+ plasticState(mappingConstitutive(2,g,i,e))%dotState(:,mappingConstitutive(1,g,i,e)) &
mySizeDotState = plasticState(mappingConstitutive(2,g,i,e))%sizeDotState
plasticState(mappingConstitutive(2,g,i,e))%state(1:mySizeDotState,mappingConstitutive(1,g,i,e)) = &
plasticState(mappingConstitutive(2,g,i,e))%subState0(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
endif
enddo; enddo; enddo
!$OMP ENDDO
@ -3187,7 +3180,7 @@ real(pReal), dimension(constitutive_maxSizeDotState) :: &
#else
mySizeDotState = plasticState(mappingConstitutive(2,g,i,e))%sizeState
mySizeDotState = plasticState(mappingConstitutive(2,g,i,e))%sizeDotState
dot_prod12 = dot_product( plasticState(mappingConstitutive(2,g,i,e))%dotState(:,mappingConstitutive(1,g,i,e)) &
- plasticState(mappingConstitutive(2,g,i,e))%previousDotState(:,mappingConstitutive(1,g,i,e)), &
plasticState(mappingConstitutive(2,g,i,e))%previousDotState(:,mappingConstitutive(1,g,i,e)) &
@ -3212,17 +3205,22 @@ real(pReal), dimension(constitutive_maxSizeDotState) :: &
endif
! --- get residui ---
mySizeDotState = plasticState(mappingConstitutive(2,g,i,e))%sizeState
stateResiduum(1:mySizeDotState) = plasticState(mappingConstitutive(2,g,i,e))%state(:,mappingConstitutive(1,g,i,e)) &
- plasticState(mappingConstitutive(2,g,i,e))%subState0(:,mappingConstitutive(1,g,i,e)) &
-(plasticState(mappingConstitutive(2,g,i,e))%dotState(:,mappingConstitutive(1,g,i,e)) * &
statedamper &
+ plasticState(mappingConstitutive(2,g,i,e))%previousDotState(:,mappingConstitutive(1,g,i,e)) &
mySizeDotState = plasticState(mappingConstitutive(2,g,i,e))%sizeDotState
stateResiduum(1:mySizeDotState) = plasticState(mappingConstitutive(2,g,i,e))% &
state(1:mySizeDotState,mappingConstitutive(1,g,i,e)) &
- plasticState(mappingConstitutive(2,g,i,e))% &
subState0(1:mySizeDotState,mappingConstitutive(1,g,i,e)) &
-(plasticState(mappingConstitutive(2,g,i,e))% &
dotState(1:mySizeDotState,mappingConstitutive(1,g,i,e)) * &
statedamper &
+ plasticState(mappingConstitutive(2,g,i,e))% &
previousDotState(1:mySizeDotState,mappingConstitutive(1,g,i,e)) &
* (1.0_pReal - statedamper)) * crystallite_subdt(g,i,e)
! --- correct state with residuum ---
tempState(1:mySizeDotState) = plasticState(mappingConstitutive(2,g,i,e))%state(:,mappingConstitutive(1,g,i,e)) &
tempState(1:mySizeDotState) = plasticState(mappingConstitutive(2,g,i,e))% &
state(1:mySizeDotState,mappingConstitutive(1,g,i,e)) &
- stateResiduum(1:mySizeDotState) ! need to copy to local variable, since we cant flush a pointer in openmp
#endif
@ -3389,13 +3387,14 @@ logical function crystallite_stateJump(g,i,e)
call constitutive_collectDeltaState(crystallite_Tstar_v(1:6,g,i,e), g,i,e)
#endif
#ifdef NEWSTATE
mySizeDotState = plasticState(mappingConstitutive(2,g,i,e))%sizeDotState
if( any(plasticState(mappingConstitutive(2,g,i,e))%deltaState(:,mappingConstitutive(1,g,i,e)) &
/= plasticState(mappingConstitutive(2,g,i,e))%deltaState(:,mappingConstitutive(1,g,i,e)))) then ! NaN occured in dotState
return
endif
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))%deltaState(:,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))%deltaState(1:mySizeDotState,mappingConstitutive(1,g,i,e))
#else
mySizeDotState = constitutive_sizeDotState(g,i,e)
if (any(constitutive_deltaState(g,i,e)%p(1:mySizeDotState) &