* 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
This commit is contained in:
Christoph Kords 2010-10-26 13:16:37 +00:00
parent d965f14f90
commit 366d52bd71
8 changed files with 186 additions and 194 deletions

View File

@ -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

View File

@ -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,7 +182,7 @@ 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
@ -205,7 +205,7 @@ 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_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))
@ -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

View File

@ -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,7 +158,7 @@ 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
@ -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

View File

@ -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

View File

@ -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

View File

@ -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

View File

@ -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,12 +1349,9 @@ 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
endif
enddo; enddo; enddo
!$OMPEND PARALLEL DO
@ -1364,10 +1361,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)
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
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
@ -1381,8 +1375,6 @@ relTemperatureResiduum = 0.0_pReal
!$OMPEND CRITICAL
endif
endif
endif
endif
enddo; enddo; enddo
!$OMPEND PARALLEL DO
@ -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,7 +1610,7 @@ relTemperatureResiduum = 0.0_pReal
! --- contribution of heun step to absolute residui ---
stateResiduum(:,g,i,e) = stateResiduum(:,g,i,e) &
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)
@ -1624,11 +1618,10 @@ relTemperatureResiduum = 0.0_pReal
! --- 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,7 +1642,9 @@ 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
@ -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

View File

@ -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 ###