From 366d52bd712b2092190dcf70972bc21e4592b9d4 Mon Sep 17 00:00:00 2001 From: Christoph Kords Date: Tue, 26 Oct 2010 13:16:37 +0000 Subject: [PATCH] * former "relevantRho" and "relevantResistance" is renamed to "atol_rho" and "atol_resistance" and now used as an absolute tolerance for the state residuum. before it was used rather as significant state, so whenever the state dropped below that value it was considered converged. (In dislotwin and titanmod constitutive law there is only one value atol_rho which is used for all the states, though the state array consists not only of densities. We need further parameters here!) * somewhat simplified convergence check for adaptive euler and runge-kutta --- code/constitutive.f90 | 33 +++++----- code/constitutive_dislotwin.f90 | 74 ++++++++++----------- code/constitutive_j2.f90 | 72 ++++++++++----------- code/constitutive_nonlocal.f90 | 46 +++++++------- code/constitutive_phenopowerlaw.f90 | 24 +++---- code/constitutive_titanmod.f90 | 20 +++--- code/crystallite.f90 | 99 ++++++++++++++--------------- code/material.config | 12 ++-- 8 files changed, 186 insertions(+), 194 deletions(-) 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 ###