diff --git a/code/constitutive.f90 b/code/constitutive.f90 index b685e0639..f78494fc6 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -23,7 +23,7 @@ MODULE constitutive constitutive_previousDotState2,&! pointer array to 2nd previous evolution of current microstructure constitutive_dotState_backup, & ! pointer array to backed up evolution of current microstructure constitutive_RK4dotState, & ! pointer array to evolution of microstructure defined by classical Runge-Kutta method - constitutive_relevantState ! relevant state values + constitutive_aTolState ! pointer array to absolute state tolerance type(p_vec), dimension(:,:,:,:), allocatable :: constitutive_RKCK45dotState ! pointer array to evolution of microstructure used by Cash-Karp Runge-Kutta method integer(pInt), dimension(:,:,:), allocatable :: constitutive_sizeDotState, & ! size of dotState array constitutive_sizeState, & ! size of state array per grain @@ -125,7 +125,7 @@ subroutine constitutive_init() allocate(constitutive_state_backup(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) allocate(constitutive_dotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) allocate(constitutive_dotState_backup(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) - allocate(constitutive_relevantState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) + allocate(constitutive_aTolState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) allocate(constitutive_sizeDotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_sizeDotState = 0_pInt allocate(constitutive_sizeState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_sizeState = 0_pInt allocate(constitutive_sizePostResults(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)); constitutive_sizePostResults = 0_pInt @@ -152,7 +152,7 @@ subroutine constitutive_init() allocate(constitutive_subState0(g,i,e)%p(constitutive_j2_sizeState(myInstance))) allocate(constitutive_state(g,i,e)%p(constitutive_j2_sizeState(myInstance))) allocate(constitutive_state_backup(g,i,e)%p(constitutive_j2_sizeState(myInstance))) - allocate(constitutive_relevantState(g,i,e)%p(constitutive_j2_sizeState(myInstance))) + allocate(constitutive_aTolState(g,i,e)%p(constitutive_j2_sizeState(myInstance))) allocate(constitutive_dotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) if (integrator == 1 .or. integratorStiffness == 1) then @@ -167,7 +167,7 @@ subroutine constitutive_init() enddo endif constitutive_state0(g,i,e)%p = constitutive_j2_stateInit(myInstance) - constitutive_relevantState(g,i,e)%p = constitutive_j2_relevantState(myInstance) + constitutive_aTolState(g,i,e)%p = constitutive_j2_aTolState(myInstance) constitutive_sizeState(g,i,e) = constitutive_j2_sizeState(myInstance) constitutive_sizeDotState(g,i,e) = constitutive_j2_sizeDotState(myInstance) constitutive_sizePostResults(g,i,e) = constitutive_j2_sizePostResults(myInstance) @@ -178,7 +178,7 @@ subroutine constitutive_init() allocate(constitutive_subState0(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) allocate(constitutive_state(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) allocate(constitutive_state_backup(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) - allocate(constitutive_relevantState(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) + allocate(constitutive_aTolState(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) allocate(constitutive_dotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) if (integrator == 1 .or. integratorStiffness == 1) then @@ -193,7 +193,7 @@ subroutine constitutive_init() enddo endif constitutive_state0(g,i,e)%p = constitutive_phenopowerlaw_stateInit(myInstance) - constitutive_relevantState(g,i,e)%p = constitutive_phenopowerlaw_relevantState(myInstance) + constitutive_aTolState(g,i,e)%p = constitutive_phenopowerlaw_aTolState(myInstance) constitutive_sizeState(g,i,e) = constitutive_phenopowerlaw_sizeState(myInstance) constitutive_sizeDotState(g,i,e) = constitutive_phenopowerlaw_sizeDotState(myInstance) constitutive_sizePostResults(g,i,e) = constitutive_phenopowerlaw_sizePostResults(myInstance) @@ -204,7 +204,7 @@ subroutine constitutive_init() allocate(constitutive_subState0(g,i,e)%p(constitutive_titanmod_sizeState(myInstance))) allocate(constitutive_state(g,i,e)%p(constitutive_titanmod_sizeState(myInstance))) allocate(constitutive_state_backup(g,i,e)%p(constitutive_titanmod_sizeState(myInstance))) - allocate(constitutive_relevantState(g,i,e)%p(constitutive_titanmod_sizeState(myInstance))) + allocate(constitutive_aTolState(g,i,e)%p(constitutive_titanmod_sizeState(myInstance))) allocate(constitutive_dotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) if (integrator == 1 .or. integratorStiffness == 1) then @@ -219,7 +219,7 @@ subroutine constitutive_init() enddo endif constitutive_state0(g,i,e)%p = constitutive_titanmod_stateInit(myInstance) - constitutive_relevantState(g,i,e)%p = constitutive_titanmod_relevantState(myInstance) + constitutive_aTolState(g,i,e)%p = constitutive_titanmod_aTolState(myInstance) constitutive_sizeState(g,i,e) = constitutive_titanmod_sizeState(myInstance) constitutive_sizeDotState(g,i,e) = constitutive_titanmod_sizeDotState(myInstance) constitutive_sizePostResults(g,i,e) = constitutive_titanmod_sizePostResults(myInstance) @@ -230,7 +230,7 @@ subroutine constitutive_init() allocate(constitutive_subState0(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance))) allocate(constitutive_state(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance))) allocate(constitutive_state_backup(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance))) - allocate(constitutive_relevantState(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance))) + allocate(constitutive_aTolState(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance))) allocate(constitutive_dotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) if (integrator == 1 .or. integratorStiffness == 1) then @@ -245,7 +245,7 @@ subroutine constitutive_init() enddo endif constitutive_state0(g,i,e)%p = constitutive_dislotwin_stateInit(myInstance) - constitutive_relevantState(g,i,e)%p = constitutive_dislotwin_relevantState(myInstance) + constitutive_aTolState(g,i,e)%p = constitutive_dislotwin_aTolState(myInstance) constitutive_sizeState(g,i,e) = constitutive_dislotwin_sizeState(myInstance) constitutive_sizeDotState(g,i,e) = constitutive_dislotwin_sizeDotState(myInstance) constitutive_sizePostResults(g,i,e) = constitutive_dislotwin_sizePostResults(myInstance) @@ -256,7 +256,7 @@ subroutine constitutive_init() allocate(constitutive_subState0(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance))) allocate(constitutive_state(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance))) allocate(constitutive_state_backup(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance))) - allocate(constitutive_relevantState(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance))) + allocate(constitutive_aTolState(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance))) allocate(constitutive_dotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) if (integrator == 1 .or. integratorStiffness == 1) then @@ -271,7 +271,7 @@ subroutine constitutive_init() enddo endif constitutive_state0(g,i,e)%p = constitutive_nonlocal_stateInit(myInstance) - constitutive_relevantState(g,i,e)%p = constitutive_nonlocal_relevantState(myInstance) + constitutive_aTolState(g,i,e)%p = constitutive_nonlocal_aTolState(myInstance) constitutive_sizeState(g,i,e) = constitutive_nonlocal_sizeState(myInstance) constitutive_sizeDotState(g,i,e) = constitutive_nonlocal_sizeDotState(myInstance) constitutive_sizePostResults(g,i,e) = constitutive_nonlocal_sizePostResults(myInstance) @@ -298,7 +298,7 @@ subroutine constitutive_init() write(6,'(a32,x,7(i5,x))') 'constitutive_partionedState0: ', shape(constitutive_partionedState0) write(6,'(a32,x,7(i5,x))') 'constitutive_subState0: ', shape(constitutive_subState0) write(6,'(a32,x,7(i5,x))') 'constitutive_state: ', shape(constitutive_state) - write(6,'(a32,x,7(i5,x))') 'constitutive_relevantState: ', shape(constitutive_relevantState) + write(6,'(a32,x,7(i5,x))') 'constitutive_aTolState: ', shape(constitutive_aTolState) write(6,'(a32,x,7(i5,x))') 'constitutive_dotState: ', shape(constitutive_dotState) write(6,'(a32,x,7(i5,x))') 'constitutive_sizeState: ', shape(constitutive_sizeState) write(6,'(a32,x,7(i5,x))') 'constitutive_sizeDotState: ', shape(constitutive_sizeDotState) @@ -499,8 +499,7 @@ 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) case (constitutive_nonlocal_label) - call constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, Temperature, constitutive_state, & - constitutive_relevantState, ipc, ip, el) + call constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, Temperature, constitutive_state, ipc, ip, el) end select @@ -579,7 +578,7 @@ select case (phase_constitution(material_phase(ipc,ip,el))) case (constitutive_nonlocal_label) call constitutive_nonlocal_dotState(constitutive_dotState, Tstar_v, subTstar0_v, Fe, Fp, Temperature, subdt, & - constitutive_state, constitutive_subState0, constitutive_relevantState, subdt, & + constitutive_state, constitutive_subState0, constitutive_aTolState, subdt, & orientation, ipc, ip, el) end select @@ -712,7 +711,7 @@ function constitutive_postResults(Tstar_v, subTstar0_v, Fe, Fp, Temperature, mis case (constitutive_nonlocal_label) constitutive_postResults = constitutive_nonlocal_postResults(Tstar_v, subTstar0_v, Fe, Fp, Temperature, misorientation, & dt, subdt, constitutive_state, constitutive_subState0, & - constitutive_relevantState, constitutive_dotstate, ipc, ip, el) + constitutive_dotstate, ipc, ip, el) end select return diff --git a/code/constitutive_dislotwin.f90 b/code/constitutive_dislotwin.f90 index a704a4e60..9702aacfa 100644 --- a/code/constitutive_dislotwin.f90 +++ b/code/constitutive_dislotwin.f90 @@ -60,7 +60,7 @@ real(pReal), dimension(:), allocatable :: constitutive_dislotwin constitutive_dislotwin_Cthresholdtwin, & ! constitutive_dislotwin_SolidSolutionStrength, & ! Strength due to elements in solid solution constitutive_dislotwin_L0, & ! Length of twin nuclei in Burgers vectors - constitutive_dislotwin_relevantRho ! dislocation density considered relevant + constitutive_dislotwin_aTolRho ! absolute tolerance for integration of dislocation density real(pReal), dimension(:,:,:), allocatable :: constitutive_dislotwin_Cslip_66 ! elasticity matrix in Mandel notation for each instance real(pReal), dimension(:,:,:,:), allocatable :: constitutive_dislotwin_Ctwin_66 ! twin elasticity matrix in Mandel notation for each instance real(pReal), dimension(:,:,:,:,:), allocatable :: constitutive_dislotwin_Cslip_3333 ! elasticity matrix for each instance @@ -182,32 +182,32 @@ allocate(constitutive_dislotwin_Cmfptwin(maxNinstance)) allocate(constitutive_dislotwin_Cthresholdtwin(maxNinstance)) allocate(constitutive_dislotwin_SolidSolutionStrength(maxNinstance)) allocate(constitutive_dislotwin_L0(maxNinstance)) -allocate(constitutive_dislotwin_relevantRho(maxNinstance)) +allocate(constitutive_dislotwin_aTolRho(maxNinstance)) allocate(constitutive_dislotwin_Cslip_66(6,6,maxNinstance)) allocate(constitutive_dislotwin_Cslip_3333(3,3,3,3,maxNinstance)) -constitutive_dislotwin_CoverA = 0.0_pReal -constitutive_dislotwin_C11 = 0.0_pReal -constitutive_dislotwin_C12 = 0.0_pReal -constitutive_dislotwin_C13 = 0.0_pReal -constitutive_dislotwin_C33 = 0.0_pReal -constitutive_dislotwin_C44 = 0.0_pReal -constitutive_dislotwin_Gmod = 0.0_pReal -constitutive_dislotwin_CAtomicVolume = 0.0_pReal -constitutive_dislotwin_D0 = 0.0_pReal -constitutive_dislotwin_Qsd = 0.0_pReal -constitutive_dislotwin_GrainSize = 0.0_pReal -constitutive_dislotwin_p = 0.0_pReal -constitutive_dislotwin_q = 0.0_pReal -constitutive_dislotwin_MaxTwinFraction = 0.0_pReal -constitutive_dislotwin_r = 0.0_pReal -constitutive_dislotwin_CEdgeDipMinDistance = 0.0_pReal -constitutive_dislotwin_Cmfptwin = 0.0_pReal -constitutive_dislotwin_Cthresholdtwin = 0.0_pReal -constitutive_dislotwin_SolidSolutionStrength = 0.0_pReal -constitutive_dislotwin_L0 = 0.0_pReal -constitutive_dislotwin_relevantRho = 0.0_pReal -constitutive_dislotwin_Cslip_66 = 0.0_pReal -constitutive_dislotwin_Cslip_3333 = 0.0_pReal +constitutive_dislotwin_CoverA = 0.0_pReal +constitutive_dislotwin_C11 = 0.0_pReal +constitutive_dislotwin_C12 = 0.0_pReal +constitutive_dislotwin_C13 = 0.0_pReal +constitutive_dislotwin_C33 = 0.0_pReal +constitutive_dislotwin_C44 = 0.0_pReal +constitutive_dislotwin_Gmod = 0.0_pReal +constitutive_dislotwin_CAtomicVolume = 0.0_pReal +constitutive_dislotwin_D0 = 0.0_pReal +constitutive_dislotwin_Qsd = 0.0_pReal +constitutive_dislotwin_GrainSize = 0.0_pReal +constitutive_dislotwin_p = 0.0_pReal +constitutive_dislotwin_q = 0.0_pReal +constitutive_dislotwin_MaxTwinFraction = 0.0_pReal +constitutive_dislotwin_r = 0.0_pReal +constitutive_dislotwin_CEdgeDipMinDistance = 0.0_pReal +constitutive_dislotwin_Cmfptwin = 0.0_pReal +constitutive_dislotwin_Cthresholdtwin = 0.0_pReal +constitutive_dislotwin_SolidSolutionStrength= 0.0_pReal +constitutive_dislotwin_L0 = 0.0_pReal +constitutive_dislotwin_aTolRho = 0.0_pReal +constitutive_dislotwin_Cslip_66 = 0.0_pReal +constitutive_dislotwin_Cslip_3333 = 0.0_pReal allocate(constitutive_dislotwin_rhoEdge0(lattice_maxNslipFamily,maxNinstance)) allocate(constitutive_dislotwin_rhoEdgeDip0(lattice_maxNslipFamily,maxNinstance)) allocate(constitutive_dislotwin_burgersPerSlipFamily(lattice_maxNslipFamily,maxNinstance)) @@ -230,10 +230,10 @@ allocate(constitutive_dislotwin_interactionSlipSlip(lattice_maxNinteraction,maxN allocate(constitutive_dislotwin_interactionSlipTwin(lattice_maxNinteraction,maxNinstance)) allocate(constitutive_dislotwin_interactionTwinSlip(lattice_maxNinteraction,maxNinstance)) allocate(constitutive_dislotwin_interactionTwinTwin(lattice_maxNinteraction,maxNinstance)) -constitutive_dislotwin_interactionSlipSlip = 0.0_pReal -constitutive_dislotwin_interactionSlipTwin = 0.0_pReal -constitutive_dislotwin_interactionTwinSlip = 0.0_pReal -constitutive_dislotwin_interactionTwinTwin = 0.0_pReal +constitutive_dislotwin_interactionSlipSlip = 0.0_pReal +constitutive_dislotwin_interactionSlipTwin = 0.0_pReal +constitutive_dislotwin_interactionTwinSlip = 0.0_pReal +constitutive_dislotwin_interactionTwinTwin = 0.0_pReal !* Readout data from material.config file rewind(file) @@ -321,15 +321,15 @@ do ! read thru sections of constitutive_dislotwin_D0(i) = IO_floatValue(line,positions,2) case ('qsd') constitutive_dislotwin_Qsd(i) = IO_floatValue(line,positions,2) - case ('relevantrho') - constitutive_dislotwin_relevantRho(i) = IO_floatValue(line,positions,2) + case ('atol_rho') + constitutive_dislotwin_aTolRho(i) = IO_floatValue(line,positions,2) case ('cmfptwin') constitutive_dislotwin_Cmfptwin(i) = IO_floatValue(line,positions,2) case ('cthresholdtwin') constitutive_dislotwin_Cthresholdtwin(i) = IO_floatValue(line,positions,2) case ('solidsolutionstrength') constitutive_dislotwin_SolidSolutionStrength(i) = IO_floatValue(line,positions,2) - case ('L0') + case ('l0') constitutive_dislotwin_L0(i) = IO_floatValue(line,positions,2) case ('cedgedipmindistance') constitutive_dislotwin_CEdgeDipMinDistance(i) = IO_floatValue(line,positions,2) @@ -378,7 +378,7 @@ enddo if (constitutive_dislotwin_CAtomicVolume(i) <= 0.0_pReal) call IO_error(230) if (constitutive_dislotwin_D0(i) <= 0.0_pReal) call IO_error(231) if (constitutive_dislotwin_Qsd(i) <= 0.0_pReal) call IO_error(232) - if (constitutive_dislotwin_relevantRho(i) <= 0.0_pReal) call IO_error(233) + if (constitutive_dislotwin_aTolRho(i) <= 0.0_pReal) call IO_error(233) !* Determine total number of active slip or twin systems constitutive_dislotwin_Nslip(:,i) = min(lattice_NslipSystem(:,myStructure),constitutive_dislotwin_Nslip(:,i)) @@ -670,18 +670,18 @@ return end function -pure function constitutive_dislotwin_relevantState(myInstance) +pure function constitutive_dislotwin_aTolState(myInstance) !********************************************************************* -!* relevant microstructural state * +!* absolute state tolerance * !********************************************************************* use prec, only: pReal, pInt implicit none !* Input-Output variables integer(pInt), intent(in) :: myInstance -real(pReal), dimension(constitutive_dislotwin_sizeState(myInstance)) :: constitutive_dislotwin_relevantState +real(pReal), dimension(constitutive_dislotwin_sizeState(myInstance)) :: constitutive_dislotwin_aTolState -constitutive_dislotwin_relevantState = constitutive_dislotwin_relevantRho(myInstance) +constitutive_dislotwin_aTolState = constitutive_dislotwin_aTolRho(myInstance) return endfunction diff --git a/code/constitutive_j2.f90 b/code/constitutive_j2.f90 index 27044440f..5e978c658 100644 --- a/code/constitutive_j2.f90 +++ b/code/constitutive_j2.f90 @@ -1,25 +1,25 @@ !* $Id$ !***************************************************** -!* Module: CONSTITUTIVE_J2 * +!* Module: CONSTITUTIVE_J2 * !***************************************************** !* contains: * !* - constitutive equations * !* - parameters definition * !***************************************************** -! [Alu] -! constitution j2 -! (output) flowstress -! (output) strainrate -! c11 110.9e9 # (3 C11 + 2 C12 + 2 C44) / 5 ... with C44 = C11-C12 !! -! c12 58.34e9 # (1 C11 + 4 C12 - 1 C44) / 5 -! taylorfactor 3 -! tau0 31e6 -! gdot0 0.001 -! n 20 -! h0 75e6 -! tausat 63e6 -! w0 2.25 +! [Alu] +! constitution j2 +! (output) flowstress +! (output) strainrate +! c11 110.9e9 # (3 C11 + 2 C12 + 2 C44) / 5 ... with C44 = C11-C12 !! +! c12 58.34e9 # (1 C11 + 4 C12 - 1 C44) / 5 +! taylorfactor 3 +! tau0 31e6 +! gdot0 0.001 +! n 20 +! h0 75e6 +! tausat 63e6 +! w0 2.25 MODULE constitutive_j2 @@ -45,7 +45,7 @@ MODULE constitutive_j2 real(pReal), dimension(:), allocatable :: constitutive_j2_h0 real(pReal), dimension(:), allocatable :: constitutive_j2_tausat real(pReal), dimension(:), allocatable :: constitutive_j2_w0 - real(pReal), dimension(:), allocatable :: constitutive_j2_relevantResistance + real(pReal), dimension(:), allocatable :: constitutive_j2_aTolResistance CONTAINS @@ -101,7 +101,7 @@ subroutine constitutive_j2_init(file) allocate(constitutive_j2_h0(maxNinstance)) ; constitutive_j2_h0 = 0.0_pReal allocate(constitutive_j2_tausat(maxNinstance)) ; constitutive_j2_tausat = 0.0_pReal allocate(constitutive_j2_w0(maxNinstance)) ; constitutive_j2_w0 = 0.0_pReal - allocate(constitutive_j2_relevantResistance(maxNinstance)) ; constitutive_j2_relevantResistance = 0.0_pReal + allocate(constitutive_j2_aTolResistance(maxNinstance)) ; constitutive_j2_aTolResistance = 0.0_pReal rewind(file) line = '' @@ -145,8 +145,8 @@ subroutine constitutive_j2_init(file) constitutive_j2_w0(i) = IO_floatValue(line,positions,2) case ('taylorfactor') constitutive_j2_fTaylor(i) = IO_floatValue(line,positions,2) - case ('relevantresistance') - constitutive_j2_relevantResistance(i) = IO_floatValue(line,positions,2) + case ('atol_resistance') + constitutive_j2_aTolResistance(i) = IO_floatValue(line,positions,2) end select endif enddo @@ -158,25 +158,25 @@ subroutine constitutive_j2_init(file) if (constitutive_j2_tausat(i) <= 0.0_pReal) call IO_error(213) if (constitutive_j2_w0(i) <= 0.0_pReal) call IO_error(241) if (constitutive_j2_fTaylor(i) <= 0.0_pReal) call IO_error(240) - if (constitutive_j2_relevantResistance(i) <= 0.0_pReal) call IO_error(242) + if (constitutive_j2_aTolResistance(i) <= 0.0_pReal) call IO_error(242) enddo do i = 1,maxNinstance do j = 1,maxval(phase_Noutput) - select case(constitutive_j2_output(j,i)) - case('flowstress') - mySize = 1_pInt - case('strainrate') - mySize = 1_pInt - case default - mySize = 0_pInt - end select - - if (mySize > 0_pInt) then ! any meaningful output found - constitutive_j2_sizePostResult(j,i) = mySize - constitutive_j2_sizePostResults(i) = & - constitutive_j2_sizePostResults(i) + mySize - endif + select case(constitutive_j2_output(j,i)) + case('flowstress') + mySize = 1_pInt + case('strainrate') + mySize = 1_pInt + case default + mySize = 0_pInt + end select + + if (mySize > 0_pInt) then ! any meaningful output found + constitutive_j2_sizePostResult(j,i) = mySize + constitutive_j2_sizePostResults(i) = & + constitutive_j2_sizePostResults(i) + mySize + endif enddo constitutive_j2_sizeDotState(i) = 1 @@ -218,7 +218,7 @@ endfunction !********************************************************************* !* relevant microstructural state * !********************************************************************* -pure function constitutive_j2_relevantState(myInstance) +pure function constitutive_j2_aTolState(myInstance) use prec, only: pReal, & pInt @@ -229,11 +229,11 @@ integer(pInt), intent(in) :: myInstance ! number specifyin !*** output variables real(pReal), dimension(constitutive_j2_sizeState(myInstance)) :: & - constitutive_j2_relevantState ! relevant state values for the current instance of this constitution + constitutive_j2_aTolState ! relevant state values for the current instance of this constitution !*** local variables -constitutive_j2_relevantState = constitutive_j2_relevantResistance(myInstance) +constitutive_j2_aTolState = constitutive_j2_aTolResistance(myInstance) endfunction diff --git a/code/constitutive_nonlocal.f90 b/code/constitutive_nonlocal.f90 index 206da79d9..15a1df1fd 100644 --- a/code/constitutive_nonlocal.f90 +++ b/code/constitutive_nonlocal.f90 @@ -58,7 +58,7 @@ real(pReal), dimension(:), allocatable :: constitutive_nonlocal_ constitutive_nonlocal_atomicVolume, & ! atomic volume constitutive_nonlocal_D0, & ! prefactor for self-diffusion coefficient constitutive_nonlocal_Qsd, & ! activation enthalpy for diffusion - constitutive_nonlocal_relevantRho, & ! dislocation density considered relevant + constitutive_nonlocal_aTolRho, & ! absolute tolerance for dislocation density in state integration constitutive_nonlocal_R ! cutoff radius for dislocation stress real(pReal), dimension(:,:,:), allocatable :: constitutive_nonlocal_Cslip_66 ! elasticity matrix in Mandel notation for each instance real(pReal), dimension(:,:,:,:,:), allocatable :: constitutive_nonlocal_Cslip_3333 ! elasticity matrix for each instance @@ -213,7 +213,7 @@ allocate(constitutive_nonlocal_Q0(maxNinstance)) allocate(constitutive_nonlocal_atomicVolume(maxNinstance)) allocate(constitutive_nonlocal_D0(maxNinstance)) allocate(constitutive_nonlocal_Qsd(maxNinstance)) -allocate(constitutive_nonlocal_relevantRho(maxNinstance)) +allocate(constitutive_nonlocal_aTolRho(maxNinstance)) allocate(constitutive_nonlocal_Cslip_66(6,6,maxNinstance)) allocate(constitutive_nonlocal_Cslip_3333(3,3,3,3,maxNinstance)) allocate(constitutive_nonlocal_R(maxNinstance)) @@ -228,7 +228,7 @@ constitutive_nonlocal_Q0 = 0.0_pReal constitutive_nonlocal_atomicVolume = 0.0_pReal constitutive_nonlocal_D0 = 0.0_pReal constitutive_nonlocal_Qsd = 0.0_pReal -constitutive_nonlocal_relevantRho = 0.0_pReal +constitutive_nonlocal_aTolRho = 0.0_pReal constitutive_nonlocal_nu = 0.0_pReal constitutive_nonlocal_Cslip_66 = 0.0_pReal constitutive_nonlocal_Cslip_3333 = 0.0_pReal @@ -335,8 +335,8 @@ do constitutive_nonlocal_D0(i) = IO_floatValue(line,positions,2) case('qsd') constitutive_nonlocal_Qsd(i) = IO_floatValue(line,positions,2) - case('relevantrho') - constitutive_nonlocal_relevantRho(i) = IO_floatValue(line,positions,2) + case('atol_rho') + constitutive_nonlocal_aTolRho(i) = IO_floatValue(line,positions,2) case ('interaction_slipslip') forall (it = 1:lattice_maxNinteraction) constitutive_nonlocal_interactionSlipSlip(it,i) = IO_floatValue(line,positions,1+it) end select @@ -379,7 +379,7 @@ enddo if (constitutive_nonlocal_atomicVolume(i) <= 0.0_pReal) call IO_error(230) if (constitutive_nonlocal_D0(i) <= 0.0_pReal) call IO_error(231) if (constitutive_nonlocal_Qsd(i) <= 0.0_pReal) call IO_error(232) - if (constitutive_nonlocal_relevantRho(i) <= 0.0_pReal) call IO_error(233) + if (constitutive_nonlocal_aTolRho(i) <= 0.0_pReal) call IO_error(233) !*** determine total number of active slip systems @@ -686,9 +686,9 @@ endfunction !********************************************************************* -!* relevant microstructural state * +!* absolute state tolerance * !********************************************************************* -pure function constitutive_nonlocal_relevantState(myInstance) +pure function constitutive_nonlocal_aTolState(myInstance) use prec, only: pReal, & pInt @@ -699,11 +699,11 @@ integer(pInt), intent(in) :: myInstance ! number specifyin !*** output variables real(pReal), dimension(constitutive_nonlocal_sizeState(myInstance)) :: & - constitutive_nonlocal_relevantState ! relevant state values for the current instance of this constitution + constitutive_nonlocal_aTolState ! absolute state tolerance for the current instance of this constitution !*** local variables -constitutive_nonlocal_relevantState = constitutive_nonlocal_relevantRho(myInstance) +constitutive_nonlocal_aTolState = constitutive_nonlocal_aTolRho(myInstance) endfunction @@ -1097,7 +1097,7 @@ endsubroutine !********************************************************************* !* calculates plastic velocity gradient and its tangent * !********************************************************************* -subroutine constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar99, Tstar_v, Temperature, state, relevantState, g, ip, el) +subroutine constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar99, Tstar_v, Temperature, state, g, ip, el) use prec, only: pReal, & pInt, & @@ -1124,8 +1124,7 @@ integer(pInt), intent(in) :: g, & ! curren el ! current element number real(pReal), intent(in) :: Temperature ! temperature type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & - state, & ! microstructural state - relevantState ! relevant microstructural state + state ! microstructural state real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola-Kirchhoff stress in Mandel notation !*** output variables @@ -1180,8 +1179,8 @@ call constitutive_nonlocal_kinetics(Tstar_v, Temperature, state, g, ip, el, dv_d !*** Calculation of gdot and its tangent -forall (s = 1:ns, t = 1:4, rhoSgl(s,t) > relevantState(g,ip,el)%p((t-1)*ns+s)) & ! no shear rate contribution for densities below relevant state - gdot(s,t) = rhoSgl(s,t) * constitutive_nonlocal_burgersPerSlipSystem(s,myInstance) * constitutive_nonlocal_v(s,t,g,ip,el) +forall (t = 1:4) & + gdot(:,t) = rhoSgl(:,t) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) * constitutive_nonlocal_v(:,t,g,ip,el) gdotTotal = sum(gdot,2) dgdotTotal_dtau = sum(rhoSgl,2) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) * dv_dtau @@ -1220,7 +1219,7 @@ endsubroutine !* rate of change of microstructure * !********************************************************************* subroutine constitutive_nonlocal_dotState(dotState, Tstar_v, previousTstar_v, Fe, Fp, Temperature, dt_previous, & - state, previousState, relevantState, timestep, orientation, g,ip,el) + state, previousState, aTolState, timestep, orientation, g,ip,el) use prec, only: pReal, & pInt, & @@ -1282,7 +1281,7 @@ real(pReal), dimension(4,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & state, & ! current microstructural state previousState, & ! previous microstructural state - relevantState ! relevant microstructural state + aTolState ! absolute state tolerance !*** input/output variables type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(inout) :: & @@ -1708,7 +1707,7 @@ do i = 1,10*ns if (i > 4*ns .and. i <= 8*ns) & ! skip immobile densities continue if (previousState(g,ip,el)%p(i) + dotState(g,ip,el)%p(i)*timestep < 0.0_pReal) then ! if single mobile densities become negative... - if (previousState(g,ip,el)%p(i) < relevantState(g,ip,el)%p(i)) then ! ... and density is already below relevance... + if (previousState(g,ip,el)%p(i) < aTolState(g,ip,el)%p(i)) then ! ... and density is already below absolute tolerance... dotState(g,ip,el)%p(i) = - previousState(g,ip,el)%p(i) / timestep ! ... set dotState to zero else ! ... otherwise... if (verboseDebugger) then @@ -1869,8 +1868,8 @@ do n = 1,FE_NipNeighbors(mesh_element(2,e)) compatibilityMax = maxval(abs(constitutive_nonlocal_compatibility(1,:,s1,n,i,e)), compatibilityMask) compatibilityMaxCount = dble(count(abs(constitutive_nonlocal_compatibility(1,:,s1,n,i,e)) == compatibilityMax)) where (abs(constitutive_nonlocal_compatibility(1,:,s1,n,i,e)) >= compatibilityMax) compatibilityMask = .false. - if (compatibilitySum + compatibilityMax * compatibilityMaxCount > 1.0_pReal) & - where (abs(constitutive_nonlocal_compatibility(:,:,s1,n,i,e)) == compatibilityMax) & + if (compatibilitySum + compatibilityMax * compatibilityMaxCount > 1.0_pReal) & ! if compatibility sum exceeds 1... + where (abs(constitutive_nonlocal_compatibility(:,:,s1,n,i,e)) == compatibilityMax) & ! ... equally distribute what is left constitutive_nonlocal_compatibility(:,:,s1,n,i,e) = sign((1.0_pReal - compatibilitySum) / compatibilityMaxCount, & constitutive_nonlocal_compatibility(:,:,s1,n,i,e)) compatibilitySum = compatibilitySum + compatibilityMaxCount * compatibilityMax @@ -1935,7 +1934,7 @@ endfunction !* return array of constitutive results * !********************************************************************* function constitutive_nonlocal_postResults(Tstar_v, previousTstar_v, Fe, Fp, Temperature, disorientation, dt, dt_previous, & - state, previousState, relevantState, dotState, g,ip,el) + state, previousState, dotState, g,ip,el) use prec, only: pReal, & pInt, & @@ -1988,7 +1987,6 @@ real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & state, & ! current microstructural state previousState, & ! previous microstructural state - relevantState, & dotState ! evolution rate of microstructural state !*** output variables @@ -2080,8 +2078,8 @@ do t = 1,4 enddo enddo -forall (s = 1:ns, t = 1:4, rhoSgl(s,t) > relevantState(g,ip,el)%p((t-1)*ns+s)) & ! no shear rate contribution for densities below relevant state - gdot(s,t) = rhoSgl(s,t) * constitutive_nonlocal_burgersPerSlipSystem(s,myInstance) * constitutive_nonlocal_v(s,t,g,ip,el) +forall (t = 1:4) & + gdot(:,t) = rhoSgl(:,t) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) * constitutive_nonlocal_v(:,t,g,ip,el) !* calculate limits for stable dipole height and its rate of change diff --git a/code/constitutive_phenopowerlaw.f90 b/code/constitutive_phenopowerlaw.f90 index 5abf40103..79020d14e 100644 --- a/code/constitutive_phenopowerlaw.f90 +++ b/code/constitutive_phenopowerlaw.f90 @@ -111,7 +111,7 @@ MODULE constitutive_phenopowerlaw real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_w0_slip - real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_relevantResistance + real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_aTolResistance CONTAINS !**************************************** @@ -221,8 +221,8 @@ subroutine constitutive_phenopowerlaw_init(file) allocate(constitutive_phenopowerlaw_w0_slip(maxNinstance)) constitutive_phenopowerlaw_w0_slip = 0.0_pReal - allocate(constitutive_phenopowerlaw_relevantResistance(maxNinstance)) - constitutive_phenopowerlaw_relevantResistance = 0.0_pReal + allocate(constitutive_phenopowerlaw_aTolResistance(maxNinstance)) + constitutive_phenopowerlaw_aTolResistance = 0.0_pReal rewind(file) line = '' @@ -300,8 +300,8 @@ subroutine constitutive_phenopowerlaw_init(file) constitutive_phenopowerlaw_h0_twinslip(i) = IO_floatValue(line,positions,2) case ('h0_twintwin') constitutive_phenopowerlaw_h0_twintwin(i) = IO_floatValue(line,positions,2) - case ('relevantresistance') - constitutive_phenopowerlaw_relevantResistance(i) = IO_floatValue(line,positions,2) + case ('atol_resistance') + constitutive_phenopowerlaw_aTolResistance(i) = IO_floatValue(line,positions,2) case ('interaction_slipslip') forall (j = 1:lattice_maxNinteraction) & constitutive_phenopowerlaw_interaction_slipslip(j,i) = IO_floatValue(line,positions,1+j) @@ -344,8 +344,8 @@ subroutine constitutive_phenopowerlaw_init(file) any(constitutive_phenopowerlaw_Ntwin(:,i) > 0)) call IO_error(211,i) if ( constitutive_phenopowerlaw_n_twin(i) <= 0.0_pReal .and. & any(constitutive_phenopowerlaw_Ntwin(:,i) > 0)) call IO_error(212,i) - if (constitutive_phenopowerlaw_relevantResistance(i) <= 0.0_pReal) & - constitutive_phenopowerlaw_relevantResistance(i) = 1.0_pReal ! default absolute tolerance 1 Pa + if (constitutive_phenopowerlaw_aTolResistance(i) <= 0.0_pReal) & + constitutive_phenopowerlaw_aTolResistance(i) = 1.0_pReal ! default absolute tolerance 1 Pa enddo @@ -522,9 +522,9 @@ endfunction !********************************************************************* -!* relevant microstructural state * +!* absolute state tolerance * !********************************************************************* -pure function constitutive_phenopowerlaw_relevantState(myInstance) +pure function constitutive_phenopowerlaw_aTolState(myInstance) use prec, only: pReal, & pInt @@ -535,11 +535,11 @@ integer(pInt), intent(in) :: myInstance ! number specifyin !*** output variables real(pReal), dimension(constitutive_phenopowerlaw_sizeState(myInstance)) :: & - constitutive_phenopowerlaw_relevantState ! relevant state values for the current instance of this constitution + constitutive_phenopowerlaw_aTolState ! relevant state values for the current instance of this constitution !*** local variables -constitutive_phenopowerlaw_relevantState = constitutive_phenopowerlaw_relevantResistance(myInstance) +constitutive_phenopowerlaw_aTolState = constitutive_phenopowerlaw_aTolResistance(myInstance) endfunction @@ -612,7 +612,7 @@ subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temp use debug, only: debugger use math, only: math_Plain3333to99 use lattice, only: lattice_Sslip,lattice_Sslip_v,lattice_Stwin,lattice_Stwin_v, lattice_maxNslipFamily, lattice_maxNtwinFamily, & - lattice_NslipSystem,lattice_NtwinSystem + lattice_NslipSystem,lattice_NtwinSystem use mesh, only: mesh_NcpElems,mesh_maxNips use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance diff --git a/code/constitutive_titanmod.f90 b/code/constitutive_titanmod.f90 index 193ddc427..1d69535c8 100644 --- a/code/constitutive_titanmod.f90 +++ b/code/constitutive_titanmod.f90 @@ -65,7 +65,7 @@ real(pReal), dimension(:), allocatable :: constitutive_titanmod_ constitutive_titanmod_CEdgeDipMinDistance, & ! Not being used constitutive_titanmod_Cmfptwin, & ! Not being used constitutive_titanmod_Cthresholdtwin, & ! Not being used - constitutive_titanmod_relevantRho ! dislocation density considered relevant + constitutive_titanmod_aTolRho ! absolute tolerance for integration of dislocation density real(pReal), dimension(:,:,:), allocatable :: constitutive_titanmod_Cslip_66 ! elasticity matrix in Mandel notation for each instance real(pReal), dimension(:,:,:,:), allocatable :: constitutive_titanmod_Ctwin_66 ! twin elasticity matrix in Mandel notation for each instance real(pReal), dimension(:,:,:,:,:), allocatable :: constitutive_titanmod_Cslip_3333 ! elasticity matrix for each instance @@ -230,7 +230,7 @@ allocate(constitutive_titanmod_r(maxNinstance)) allocate(constitutive_titanmod_CEdgeDipMinDistance(maxNinstance)) allocate(constitutive_titanmod_Cmfptwin(maxNinstance)) allocate(constitutive_titanmod_Cthresholdtwin(maxNinstance)) -allocate(constitutive_titanmod_relevantRho(maxNinstance)) +allocate(constitutive_titanmod_aTolRho(maxNinstance)) allocate(constitutive_titanmod_Cslip_66(6,6,maxNinstance)) allocate(constitutive_titanmod_Cslip_3333(3,3,3,3,maxNinstance)) constitutive_titanmod_CoverA = 0.0_pReal @@ -251,7 +251,7 @@ constitutive_titanmod_r = 0.0_pReal constitutive_titanmod_CEdgeDipMinDistance = 0.0_pReal constitutive_titanmod_Cmfptwin = 0.0_pReal constitutive_titanmod_Cthresholdtwin = 0.0_pReal -constitutive_titanmod_relevantRho = 0.0_pReal +constitutive_titanmod_aTolRho = 0.0_pReal constitutive_titanmod_Cslip_66 = 0.0_pReal constitutive_titanmod_Cslip_3333 = 0.0_pReal allocate(constitutive_titanmod_rho_edge0(lattice_maxNslipFamily,maxNinstance)) @@ -504,8 +504,8 @@ do ! read thru sections of case ('twinhpconstant') constitutive_titanmod_twinhpconstant(i) = IO_floatValue(line,positions,2) write(6,*) tag - case ('relevantrho') - constitutive_titanmod_relevantRho(i) = IO_floatValue(line,positions,2) + case ('atol_rho') + constitutive_titanmod_aTolRho(i) = IO_floatValue(line,positions,2) write(6,*) tag case ('interactionslipslip') forall (j = 1:lattice_maxNinteraction) & @@ -577,7 +577,7 @@ write(6,*) 'Material Property reading done' ! if (any(constitutive_titanmod_interactionSlipSlip(1:maxval(lattice_interactionSlipSlip(:,:,myStructure)),i) < 1.0_pReal)) call IO_error(229) if (constitutive_titanmod_dc(i) <= 0.0_pReal) call IO_error(231) if (constitutive_titanmod_twinhpconstant(i) <= 0.0_pReal) call IO_error(232) - if (constitutive_titanmod_relevantRho(i) <= 0.0_pReal) call IO_error(233) + if (constitutive_titanmod_aTolRho(i) <= 0.0_pReal) call IO_error(233) !* Determine total number of active slip or twin systems constitutive_titanmod_Nslip(:,i) = min(lattice_NslipSystem(:,myStructure),constitutive_titanmod_Nslip(:,i)) @@ -1023,18 +1023,18 @@ constitutive_titanmod_stateInit(6*ns+nt+1:6*ns+2*nt)=resistance_twin0 return end function -pure function constitutive_titanmod_relevantState(myInstance) +pure function constitutive_titanmod_aTolState(myInstance) !********************************************************************* -!* relevant microstructural state * +!* absolute state tolerance * !********************************************************************* use prec, only: pReal, pInt implicit none !* Input-Output variables integer(pInt), intent(in) :: myInstance -real(pReal), dimension(constitutive_titanmod_sizeState(myInstance)) :: constitutive_titanmod_relevantState +real(pReal), dimension(constitutive_titanmod_sizeState(myInstance)) :: constitutive_titanmod_aTolState -constitutive_titanmod_relevantState = constitutive_titanmod_relevantRho(myInstance) +constitutive_titanmod_aTolState = constitutive_titanmod_aTolRho(myInstance) return endfunction diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 6712b788f..9a7131a1c 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -1036,7 +1036,7 @@ use material, only: homogenization_Ngrains, & use constitutive, only: constitutive_sizeDotState, & constitutive_maxSizeDotState, & constitutive_state, & - constitutive_relevantState, & + constitutive_aTolState, & constitutive_subState0, & constitutive_dotState, & constitutive_RKCK45dotState, & @@ -1078,8 +1078,6 @@ real(pReal), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: temperatureResiduum, & ! residuum from evolution in temperature relTemperatureResiduum ! relative residuum from evolution in temperature logical singleRun ! flag indicating computation for single (g,i,e) triple -logical, dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & - stateConverged ! flag indicating state convergence ! --- FILL BUTCHER TABLEAU --- @@ -1157,7 +1155,7 @@ endif !$OMP PARALLEL 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 - selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode==1) + selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, & crystallite_Fp, crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), & crystallite_orientation, g,i,e) @@ -1223,7 +1221,7 @@ do n = 1,5 !$OMP PARALLEL 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 - selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode == 1) + selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) if (crystallite_todo(g,i,e)) then if (.not. crystallite_integrateStress(mode,g,i,e,c(n))) then ! fraction of original time step if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local... @@ -1243,7 +1241,7 @@ do n = 1,5 !$OMP PARALLEL 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 - selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode==1) + selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) if (crystallite_todo(g,i,e)) then call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & @@ -1280,7 +1278,7 @@ relTemperatureResiduum = 0.0_pReal !$OMP PARALLEL 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 - selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode==1) + selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) sizeDotState = constitutive_sizeDotState(g,i,e) @@ -1309,15 +1307,17 @@ relTemperatureResiduum = 0.0_pReal + constitutive_dotState(g,i,e)%p(1:sizeDotState) * crystallite_subdt(g,i,e) crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) & + crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e) - forall (s = 1:sizeDotState, constitutive_state(g,i,e)%p(s) > constitutive_relevantState(g,i,e)%p(s)) & - relStateResiduum(s,g,i,e) = abs(stateResiduum(s,g,i,e)) / constitutive_state(g,i,e)%p(s) / rTol_crystalliteState + forall (s = 1:sizeDotState, abs(constitutive_state(g,i,e)%p(s)) > 0.0_pReal) & + relStateResiduum(s,g,i,e) = abs(stateResiduum(s,g,i,e)) / constitutive_state(g,i,e)%p(s) if (crystallite_Temperature(g,i,e) > 0) & - relTemperatureResiduum(g,i,e) = abs(temperatureResiduum(g,i,e)) / crystallite_Temperature(g,i,e) & - / rTol_crystalliteTemperature + relTemperatureResiduum(g,i,e) = abs(temperatureResiduum(g,i,e)) / crystallite_Temperature(g,i,e) ! --- state convergence --- -! stateConverged(g,i,e) = + crystallite_todo(g,i,e) = (all( relStateResiduum(:,g,i,e) < rTol_crystalliteState & + .or. abs(stateResiduum(1:sizeDotState,g,i,e)) < constitutive_aTolState(g,i,e)%p(1:sizeDotState))& + .and. relTemperatureResiduum(g,i,e) < rTol_crystalliteTemperature ) + if (verboseDebugger .and. selectiveDebugger) then !$OMP CRITICAL (write2out) @@ -1325,7 +1325,7 @@ relTemperatureResiduum = 0.0_pReal write(6,*) write(6,'(a,/,12(f12.1,x))') 'updateState: absolute residuum', stateResiduum(1:sizeDotState,g,i,e) write(6,*) - write(6,'(a,/,12(f12.1,x))') 'updateState: resid tolerance', relStateResiduum(1:sizeDotState,g,i,e) + write(6,'(a,/,12(f12.1,x))') 'updateState: resid tolerance',relStateResiduum(1:sizeDotState,g,i,e) / rTol_crystalliteState write(6,*) ! write(6,'(a)') 'updateState: RKCK45dotState' ! do j = 1,6 @@ -1349,11 +1349,8 @@ relTemperatureResiduum = 0.0_pReal !$OMP PARALLEL 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 - sizeDotState = constitutive_sizeDotState(g,i,e) - if ( all(relStateResiduum(1:sizeDotState,g,i,e) < 1.0_pReal) .and. relTemperatureResiduum(g,i,e) < 1.0_pReal ) then - call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), & - crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states - endif + call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), & + crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states endif enddo; enddo; enddo !$OMPEND PARALLEL DO @@ -1364,24 +1361,19 @@ relTemperatureResiduum = 0.0_pReal !$OMP PARALLEL 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 - selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode==1) - sizeDotState = constitutive_sizeDotState(g,i,e) - if ( all(relStateResiduum(1:sizeDotState,g,i,e) < 1.0_pReal) .and. relTemperatureResiduum(g,i,e) < 1.0_pReal ) then - - if (crystallite_integrateStress(mode,g,i,e)) then - crystallite_converged(g,i,e) = .true. ! ... converged per definitionem - crystallite_todo(g,i,e) = .false. ! ... integration done - !$OMP CRITICAL (distributionState) - debug_StateLoopDistribution(6,mode) = debug_StateLoopDistribution(6,mode) + 1 - !$OMPEND CRITICAL (distributionState) - else - if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local... - !$OMP CRITICAL - crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped - !$OMPEND CRITICAL - endif - endif - + selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) + if (crystallite_integrateStress(mode,g,i,e)) then + crystallite_converged(g,i,e) = .true. ! ... converged per definitionem + crystallite_todo(g,i,e) = .false. ! ... integration done + !$OMP CRITICAL (distributionState) + debug_StateLoopDistribution(6,mode) = debug_StateLoopDistribution(6,mode) + 1 + !$OMPEND CRITICAL (distributionState) + else + if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local... + !$OMP CRITICAL + crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped + !$OMPEND CRITICAL + endif endif endif enddo; enddo; enddo @@ -1432,7 +1424,7 @@ use material, only: homogenization_Ngrains, & use constitutive, only: constitutive_sizeDotState, & constitutive_maxSizeDotState, & constitutive_state, & - constitutive_relevantState, & + constitutive_aTolState, & constitutive_subState0, & constitutive_dotState, & constitutive_collectDotState, & @@ -1501,9 +1493,11 @@ endif ! --- DOT STATE AND TEMPERATURE (EULER INTEGRATION) --- +stateResiduum = 0.0_pReal !$OMP PARALLEL 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 + sizeDotState = constitutive_sizeDotState(g,i,e) selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode==1) call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, & @@ -1522,7 +1516,7 @@ endif crystallite_todo(g,i,e) = .false. ! ... skip this one next time endif else - stateResiduum(:,g,i,e) = - 0.5_pReal * constitutive_dotState(g,i,e)%p * crystallite_subdt(g,i,e) ! contribution to absolute residuum in state and temperature + stateResiduum(1:sizeDotState,g,i,e) = - 0.5_pReal * constitutive_dotState(g,i,e)%p * crystallite_subdt(g,i,e) ! contribution to absolute residuum in state and temperature temperatureResiduum(g,i,e) = - 0.5_pReal * crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e) endif endif @@ -1616,19 +1610,18 @@ relTemperatureResiduum = 0.0_pReal ! --- contribution of heun step to absolute residui --- - stateResiduum(:,g,i,e) = stateResiduum(:,g,i,e) & - + 0.5_pReal * constitutive_dotState(g,i,e)%p * crystallite_subdt(g,i,e) ! contribution to absolute residuum in state and temperature + stateResiduum(1:sizeDotState,g,i,e) = stateResiduum(1:sizeDotState,g,i,e) & + + 0.5_pReal * constitutive_dotState(g,i,e)%p * crystallite_subdt(g,i,e) ! contribution to absolute residuum in state and temperature temperatureResiduum(g,i,e) = temperatureResiduum(g,i,e) & - + 0.5_pReal * crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e) + + 0.5_pReal * crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e) ! --- relative residui --- - forall (s = 1:sizeDotState, constitutive_state(g,i,e)%p(s) > constitutive_relevantState(g,i,e)%p(s)) & - relStateResiduum(s,g,i,e) = abs(stateResiduum(s,g,i,e)) / constitutive_state(g,i,e)%p(s) / rTol_crystalliteState + forall (s = 1:sizeDotState, abs(constitutive_state(g,i,e)%p(s)) > 0.0_pReal) & + relStateResiduum(s,g,i,e) = abs(stateResiduum(s,g,i,e)) / constitutive_state(g,i,e)%p(s) if (crystallite_Temperature(g,i,e) > 0) & - relTemperatureResiduum(g,i,e) = abs(temperatureResiduum(g,i,e)) / crystallite_Temperature(g,i,e) & - / rTol_crystalliteTemperature + relTemperatureResiduum(g,i,e) = abs(temperatureResiduum(g,i,e)) / crystallite_Temperature(g,i,e) if (verboseDebugger .and. selectiveDebugger) then !$OMP CRITICAL (write2out) @@ -1636,7 +1629,7 @@ relTemperatureResiduum = 0.0_pReal write(6,*) write(6,'(a,/,12(f12.1,x))') 'updateState: absolute residuum', stateResiduum(1:sizeDotState,g,i,e) write(6,*) - write(6,'(a,/,12(f12.1,x))') 'updateState: resid tolerance', relStateResiduum(1:sizeDotState,g,i,e) + write(6,'(a,/,12(f12.1,x))') 'updateState: resid tolerance',relStateResiduum(1:sizeDotState,g,i,e) / rTol_crystalliteState write(6,*) write(6,'(a,/,12(e12.5,x))') 'updateState: dotState', constitutive_dotState(g,i,e)%p(1:sizeDotState) & - 2.0_pReal * stateResiduum(1:sizeDotState,g,i,e) / crystallite_subdt(g,i,e) ! calculate former dotstate from higher order solution and state residuum @@ -1649,14 +1642,16 @@ relTemperatureResiduum = 0.0_pReal ! --- converged ? --- - if ( all(relStateResiduum(1:sizeDotState,g,i,e) < 1.0_pReal) .and. relTemperatureResiduum(g,i,e) < 1.0_pReal ) then - + if ( all( relStateResiduum(:,g,i,e) < rTol_crystalliteState & + .or. abs(stateResiduum(1:sizeDotState,g,i,e)) < constitutive_aTolState(g,i,e)%p(1:sizeDotState)) & + .and. relTemperatureResiduum(g,i,e) < rTol_crystalliteTemperature ) then + crystallite_converged(g,i,e) = .true. ! ... converged per definitionem crystallite_todo(g,i,e) = .false. ! ... integration done !$OMP CRITICAL (distributionState) debug_StateLoopDistribution(6,mode) = debug_StateLoopDistribution(6,mode) + 1 !$OMPEND CRITICAL (distributionState) - + endif endif @@ -2122,7 +2117,7 @@ use constitutive, only: constitutive_dotState, & constitutive_sizeDotState, & constitutive_subState0, & constitutive_state, & - constitutive_relevantState, & + constitutive_aTolState, & constitutive_microstructure use debug, only: debugger, & selectiveDebugger, & @@ -2163,8 +2158,8 @@ endif constitutive_state(g,i,e)%p(1:mySize) = constitutive_state(g,i,e)%p(1:mySize) - residuum -! setting flag to true if state is below relative tolerance, otherwise set it to false -crystallite_updateState = all( constitutive_state(g,i,e)%p(1:mySize) < constitutive_relevantState(g,i,e)%p(1:mySize) & +! setting flag to true if residuum is below relative/absolute tolerance, otherwise set it to false +crystallite_updateState = all( abs(residuum) < constitutive_aTolState(g,i,e)%p(1:mySize) & .or. abs(residuum) < rTol_crystalliteState*abs(constitutive_state(g,i,e)%p(1:mySize)) ) if (verboseDebugger .and. selectiveDebugger) then diff --git a/code/material.config b/code/material.config index ac43c52c2..4d7a6ff05 100644 --- a/code/material.config +++ b/code/material.config @@ -94,7 +94,7 @@ n 20 h0 75e6 tausat 63e6 w0 2.25 -relevantResistance 1 +atol_resistance 1 [Aluminum_phenopowerlaw] # slip only @@ -138,7 +138,7 @@ interaction_slipslip 1 1 1.4 1.4 1.4 1.4 interaction_sliptwin 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 interaction_twinslip 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 interaction_twintwin 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 -relevantResistance 1 +atol_resistance 1 [Aluminum_nonlocal] @@ -239,7 +239,7 @@ lambda0 100 0 0 0 # prefactor for mean free path atomicVolume 1.7e-29 D0 1e-4 # prefactor for self-diffusion coefficient Qsd 2.3e-19 # activation enthalpy for seld-diffusion -relevantRho 1e3 # dislocation density considered relevant +atol_rho 1e2 # dislocation density considered relevant interaction_SlipSlip 0.122 0.122 0.625 0.07 0.137 0.122 # Dislocation interaction coefficient [BCC_Ferrite] @@ -271,7 +271,7 @@ interaction_sliptwin 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 interaction_twinslip 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 interaction_twintwin 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 w0_slip 1.0 -relevantResistance 1 +atol_resistance 1 [BCC_Martensite] constitution phenopowerlaw @@ -302,7 +302,7 @@ interaction_sliptwin 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 interaction_twinslip 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 interaction_twintwin 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 w0_slip 2.0 -relevantResistance 1 +atol_resistance 1 [TWIP steel FeMnC] @@ -343,7 +343,7 @@ pexponent 1.0 # p-exponent in glide velocity qexponent 1.0 # q-exponent in glide velocity Catomicvolume 1.0 # Adj. parameter controlling the atomic volume [in b] Cedgedipmindistance 1.0 # Adj. parameter controlling the minimum dipole distance [in b] -relevantRho 1.0e-200 +atol_rho 1.0e-200 interactionSlipSlip 0.122 0.122 0.625 0.07 0.137 0.122 # Interaction coefficients (Kubin et al. 2008) ### Twinning parameters ###