From 401b31c951a19ad20238877d83123f254dde9bc8 Mon Sep 17 00:00:00 2001 From: Luv Sharma Date: Wed, 11 Jun 2014 12:27:41 +0000 Subject: [PATCH] new state related changes --- code/constitutive.f90 | 120 ++++++---- code/constitutive_titanmod.f90 | 400 ++++++++++++++++++++++++++------- code/crystallite.f90 | 53 +++-- 3 files changed, 421 insertions(+), 152 deletions(-) diff --git a/code/constitutive.f90 b/code/constitutive.f90 index f1b868fff..8a4d4fabd 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -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))% & diff --git a/code/constitutive_titanmod.f90 b/code/constitutive_titanmod.f90 index bc9eda712..1c75cbc48 100644 --- a/code/constitutive_titanmod.f90 +++ b/code/constitutive_titanmod.f90 @@ -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 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)) diff --git a/code/crystallite.f90 b/code/crystallite.f90 index f956f52fe..68c033318 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -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) &