223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

unifying names

This commit is contained in:
Martin Diehl 2020-12-23 07:14:07 +01:00
parent 8cf1035cf3
commit fe6a82ecc1
2 changed files with 333 additions and 334 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

373
src/constitutive.f90
View File

@ -347,8 +347,8 @@ module constitutive
dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient
end subroutine constitutive_hooke_SandItsTangents
module subroutine integrateStateFPI(g,i,e)
integer, intent(in) :: e, i, g
module subroutine integrateStateFPI(co,ip,el)
integer, intent(in) :: co, ip, el
end subroutine integrateStateFPI
end interface
@ -746,19 +746,19 @@ end subroutine constitutive_allocateState
!--------------------------------------------------------------------------------------------------
!> @brief Restore data after homog cutback.
!--------------------------------------------------------------------------------------------------
subroutine constitutive_restore(i,e)
subroutine constitutive_restore(ip,el)
integer, intent(in) :: &
i, & !< integration point number
e !< element number
ip, & !< integration point number
el !< element number
integer :: &
c, & !< constituent number
co, & !< constituent number
s
do c = 1,homogenization_Nconstituents(material_homogenizationAt(e))
do s = 1, phase_Nsources(material_phaseAt(c,e))
sourceState(material_phaseAt(c,e))%p(s)%state( :,material_phasememberAt(c,i,e)) = &
sourceState(material_phaseAt(c,e))%p(s)%partitionedState0(:,material_phasememberAt(c,i,e))
do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
do s = 1, phase_Nsources(material_phaseAt(co,el))
sourceState(material_phaseAt(co,el))%p(s)%state( :,material_phasememberAt(co,ip,el)) = &
sourceState(material_phaseAt(co,el))%p(s)%partitionedState0(:,material_phasememberAt(co,ip,el))
enddo
enddo
@ -827,7 +827,6 @@ subroutine crystallite_init
iMax, & !< maximum number of integration points
eMax !< maximum number of elements
class(tNode), pointer :: &
num_crystallite, &
debug_crystallite, & ! pointer to debug options for crystallite
@ -835,6 +834,7 @@ subroutine crystallite_init
phase, &
mech
print'(/,a)', ' <<<+- crystallite init -+>>>'
debug_crystallite => config_debug%get('crystallite', defaultVal=emptyList)
@ -986,9 +986,9 @@ function crystallite_stress()
integer :: &
NiterationCrystallite, & ! number of iterations in crystallite loop
c, & !< counter in integration point component loop
i, & !< counter in integration point loop
e, & !< counter in element loop
s, p, m
ip, & !< counter in integration point loop
el, & !< counter in element loop
s, ph, me
logical, dimension(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems) :: todo !ToDo: need to set some values to false for different Ngrains
real(pReal), dimension(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems) :: subFrac !ToDo: need to set some values to false for different Ngrains
real(pReal), dimension(:,:,:,:,:), allocatable :: &
@ -1003,27 +1003,27 @@ function crystallite_stress()
!--------------------------------------------------------------------------------------------------
! initialize to starting condition
crystallite_subStep = 0.0_pReal
!$OMP PARALLEL DO PRIVATE(p,m)
elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2); do c = 1,homogenization_Nconstituents(material_homogenizationAt(e))
p = material_phaseAt(c,e)
m = material_phaseMemberAt(c,i,e)
subLi0(1:3,1:3,c,i,e) = constitutive_mech_partionedLi0(p)%data(1:3,1:3,m)
homogenizationRequestsCalculation: if (crystallite_requested(c,i,e)) then
plasticState (material_phaseAt(c,e))%subState0( :,material_phaseMemberAt(c,i,e)) = &
plasticState (material_phaseAt(c,e))%partitionedState0(:,material_phaseMemberAt(c,i,e))
!$OMP PARALLEL DO PRIVATE(ph,me)
elementLooping1: do el = FEsolving_execElem(1),FEsolving_execElem(2)
do ip = FEsolving_execIP(1),FEsolving_execIP(2); do c = 1,homogenization_Nconstituents(material_homogenizationAt(el))
ph = material_phaseAt(c,el)
me = material_phaseMemberAt(c,ip,el)
subLi0(1:3,1:3,c,ip,el) = constitutive_mech_partionedLi0(ph)%data(1:3,1:3,me)
homogenizationRequestsCalculation: if (crystallite_requested(c,ip,el)) then
plasticState (material_phaseAt(c,el))%subState0( :,material_phaseMemberAt(c,ip,el)) = &
plasticState (material_phaseAt(c,el))%partitionedState0(:,material_phaseMemberAt(c,ip,el))
do s = 1, phase_Nsources(material_phaseAt(c,e))
sourceState(material_phaseAt(c,e))%p(s)%subState0( :,material_phaseMemberAt(c,i,e)) = &
sourceState(material_phaseAt(c,e))%p(s)%partitionedState0(:,material_phaseMemberAt(c,i,e))
do s = 1, phase_Nsources(material_phaseAt(c,el))
sourceState(material_phaseAt(c,el))%p(s)%subState0( :,material_phaseMemberAt(c,ip,el)) = &
sourceState(material_phaseAt(c,el))%p(s)%partitionedState0(:,material_phaseMemberAt(c,ip,el))
enddo
crystallite_subFp0(1:3,1:3,c,i,e) = constitutive_mech_partionedFp0(p)%data(1:3,1:3,m)
crystallite_subFi0(1:3,1:3,c,i,e) = constitutive_mech_partionedFi0(p)%data(1:3,1:3,m)
crystallite_subF0(1:3,1:3,c,i,e) = crystallite_partitionedF0(1:3,1:3,c,i,e)
subFrac(c,i,e) = 0.0_pReal
crystallite_subStep(c,i,e) = 1.0_pReal/num%subStepSizeCryst
todo(c,i,e) = .true.
crystallite_converged(c,i,e) = .false. ! pretend failed step of 1/subStepSizeCryst
crystallite_subFp0(1:3,1:3,c,ip,el) = constitutive_mech_partionedFp0(ph)%data(1:3,1:3,me)
crystallite_subFi0(1:3,1:3,c,ip,el) = constitutive_mech_partionedFi0(ph)%data(1:3,1:3,me)
crystallite_subF0(1:3,1:3,c,ip,el) = crystallite_partitionedF0(1:3,1:3,c,ip,el)
subFrac(c,ip,el) = 0.0_pReal
crystallite_subStep(c,ip,el) = 1.0_pReal/num%subStepSizeCryst
todo(c,ip,el) = .true.
crystallite_converged(c,ip,el) = .false. ! pretend failed step of 1/subStepSizeCryst
endif homogenizationRequestsCalculation
enddo; enddo
enddo elementLooping1
@ -1037,71 +1037,71 @@ function crystallite_stress()
if (debugCrystallite%extensive) &
print'(a,i6)', '<< CRYST stress >> crystallite iteration ',NiterationCrystallite
#endif
!$OMP PARALLEL DO PRIVATE(formerSubStep,p,m)
elementLooping3: do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do c = 1,homogenization_Nconstituents(material_homogenizationAt(e))
p = material_phaseAt(c,e)
m = material_phaseMemberAt(c,i,e)
!$OMP PARALLEL DO PRIVATE(formerSubStep,ph,me)
elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2)
do ip = FEsolving_execIP(1),FEsolving_execIP(2)
do c = 1,homogenization_Nconstituents(material_homogenizationAt(el))
ph = material_phaseAt(c,el)
me = material_phaseMemberAt(c,ip,el)
!--------------------------------------------------------------------------------------------------
! wind forward
if (crystallite_converged(c,i,e)) then
formerSubStep = crystallite_subStep(c,i,e)
subFrac(c,i,e) = subFrac(c,i,e) + crystallite_subStep(c,i,e)
crystallite_subStep(c,i,e) = min(1.0_pReal - subFrac(c,i,e), &
num%stepIncreaseCryst * crystallite_subStep(c,i,e))
if (crystallite_converged(c,ip,el)) then
formerSubStep = crystallite_subStep(c,ip,el)
subFrac(c,ip,el) = subFrac(c,ip,el) + crystallite_subStep(c,ip,el)
crystallite_subStep(c,ip,el) = min(1.0_pReal - subFrac(c,ip,el), &
num%stepIncreaseCryst * crystallite_subStep(c,ip,el))
todo(c,i,e) = crystallite_subStep(c,i,e) > 0.0_pReal ! still time left to integrate on?
todo(c,ip,el) = crystallite_subStep(c,ip,el) > 0.0_pReal ! still time left to integrate on?
if (todo(c,i,e)) then
crystallite_subF0 (1:3,1:3,c,i,e) = crystallite_subF(1:3,1:3,c,i,e)
subLp0(1:3,1:3,c,i,e) = crystallite_Lp (1:3,1:3,c,i,e)
subLi0(1:3,1:3,c,i,e) = constitutive_mech_Li(p)%data(1:3,1:3,m)
crystallite_subFp0(1:3,1:3,c,i,e) = constitutive_mech_Fp(p)%data(1:3,1:3,m)
crystallite_subFi0(1:3,1:3,c,i,e) = constitutive_mech_Fi(p)%data(1:3,1:3,m)
plasticState( material_phaseAt(c,e))%subState0(:,material_phaseMemberAt(c,i,e)) &
= plasticState(material_phaseAt(c,e))%state( :,material_phaseMemberAt(c,i,e))
do s = 1, phase_Nsources(material_phaseAt(c,e))
sourceState( material_phaseAt(c,e))%p(s)%subState0(:,material_phaseMemberAt(c,i,e)) &
= sourceState(material_phaseAt(c,e))%p(s)%state( :,material_phaseMemberAt(c,i,e))
if (todo(c,ip,el)) then
crystallite_subF0 (1:3,1:3,c,ip,el) = crystallite_subF(1:3,1:3,c,ip,el)
subLp0(1:3,1:3,c,ip,el) = crystallite_Lp (1:3,1:3,c,ip,el)
subLi0(1:3,1:3,c,ip,el) = constitutive_mech_Li(ph)%data(1:3,1:3,me)
crystallite_subFp0(1:3,1:3,c,ip,el) = constitutive_mech_Fp(ph)%data(1:3,1:3,me)
crystallite_subFi0(1:3,1:3,c,ip,el) = constitutive_mech_Fi(ph)%data(1:3,1:3,me)
plasticState( material_phaseAt(c,el))%subState0(:,material_phaseMemberAt(c,ip,el)) &
= plasticState(material_phaseAt(c,el))%state( :,material_phaseMemberAt(c,ip,el))
do s = 1, phase_Nsources(material_phaseAt(c,el))
sourceState( material_phaseAt(c,el))%p(s)%subState0(:,material_phaseMemberAt(c,ip,el)) &
= sourceState(material_phaseAt(c,el))%p(s)%state( :,material_phaseMemberAt(c,ip,el))
enddo
endif
!--------------------------------------------------------------------------------------------------
! cut back (reduced time and restore)
else
crystallite_subStep(c,i,e) = num%subStepSizeCryst * crystallite_subStep(c,i,e)
constitutive_mech_Fp(p)%data(1:3,1:3,m) = crystallite_subFp0(1:3,1:3,c,i,e)
constitutive_mech_Fi(p)%data(1:3,1:3,m) = crystallite_subFi0(1:3,1:3,c,i,e)
crystallite_S (1:3,1:3,c,i,e) = crystallite_S0 (1:3,1:3,c,i,e)
if (crystallite_subStep(c,i,e) < 1.0_pReal) then ! actual (not initial) cutback
crystallite_Lp (1:3,1:3,c,i,e) = subLp0(1:3,1:3,c,i,e)
constitutive_mech_Li(p)%data(1:3,1:3,m) = subLi0(1:3,1:3,c,i,e)
crystallite_subStep(c,ip,el) = num%subStepSizeCryst * crystallite_subStep(c,ip,el)
constitutive_mech_Fp(ph)%data(1:3,1:3,me) = crystallite_subFp0(1:3,1:3,c,ip,el)
constitutive_mech_Fi(ph)%data(1:3,1:3,me) = crystallite_subFi0(1:3,1:3,c,ip,el)
crystallite_S (1:3,1:3,c,ip,el) = crystallite_S0 (1:3,1:3,c,ip,el)
if (crystallite_subStep(c,ip,el) < 1.0_pReal) then ! actual (not initial) cutback
crystallite_Lp (1:3,1:3,c,ip,el) = subLp0(1:3,1:3,c,ip,el)
constitutive_mech_Li(ph)%data(1:3,1:3,me) = subLi0(1:3,1:3,c,ip,el)
endif
plasticState (material_phaseAt(c,e))%state( :,material_phaseMemberAt(c,i,e)) &
= plasticState(material_phaseAt(c,e))%subState0(:,material_phaseMemberAt(c,i,e))
do s = 1, phase_Nsources(material_phaseAt(c,e))
sourceState( material_phaseAt(c,e))%p(s)%state( :,material_phaseMemberAt(c,i,e)) &
= sourceState(material_phaseAt(c,e))%p(s)%subState0(:,material_phaseMemberAt(c,i,e))
plasticState (material_phaseAt(c,el))%state( :,material_phaseMemberAt(c,ip,el)) &
= plasticState(material_phaseAt(c,el))%subState0(:,material_phaseMemberAt(c,ip,el))
do s = 1, phase_Nsources(material_phaseAt(c,el))
sourceState( material_phaseAt(c,el))%p(s)%state( :,material_phaseMemberAt(c,ip,el)) &
= sourceState(material_phaseAt(c,el))%p(s)%subState0(:,material_phaseMemberAt(c,ip,el))
enddo
! cant restore dotState here, since not yet calculated in first cutback after initialization
todo(c,i,e) = crystallite_subStep(c,i,e) > num%subStepMinCryst ! still on track or already done (beyond repair)
todo(c,ip,el) = crystallite_subStep(c,ip,el) > num%subStepMinCryst ! still on track or already done (beyond repair)
endif
!--------------------------------------------------------------------------------------------------
! prepare for integration
if (todo(c,i,e)) then
crystallite_subF(1:3,1:3,c,i,e) = crystallite_subF0(1:3,1:3,c,i,e) &
+ crystallite_subStep(c,i,e) *( crystallite_partitionedF (1:3,1:3,c,i,e) &
-crystallite_partitionedF0(1:3,1:3,c,i,e))
crystallite_Fe(1:3,1:3,c,i,e) = matmul(crystallite_subF(1:3,1:3,c,i,e), &
math_inv33(matmul(constitutive_mech_Fi(p)%data(1:3,1:3,m), &
constitutive_mech_Fp(p)%data(1:3,1:3,m))))
crystallite_subdt(c,i,e) = crystallite_subStep(c,i,e) * crystallite_dt(c,i,e)
crystallite_converged(c,i,e) = .false.
call integrateState(c,i,e)
call integrateSourceState(c,i,e)
if (todo(c,ip,el)) then
crystallite_subF(1:3,1:3,c,ip,el) = crystallite_subF0(1:3,1:3,c,ip,el) &
+ crystallite_subStep(c,ip,el) *( crystallite_partitionedF (1:3,1:3,c,ip,el) &
-crystallite_partitionedF0(1:3,1:3,c,ip,el))
crystallite_Fe(1:3,1:3,c,ip,el) = matmul(crystallite_subF(1:3,1:3,c,ip,el), &
math_inv33(matmul(constitutive_mech_Fi(ph)%data(1:3,1:3,me), &
constitutive_mech_Fp(ph)%data(1:3,1:3,me))))
crystallite_subdt(c,ip,el) = crystallite_subStep(c,ip,el) * crystallite_dt(c,ip,el)
crystallite_converged(c,ip,el) = .false.
call integrateState(c,ip,el)
call integrateSourceState(c,ip,el)
endif
enddo
@ -1117,9 +1117,9 @@ function crystallite_stress()
! return whether converged or not
crystallite_stress = .false.
elementLooping5: do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
crystallite_stress(i,e) = all(crystallite_converged(:,i,e))
elementLooping5: do el = FEsolving_execElem(1),FEsolving_execElem(2)
do ip = FEsolving_execIP(1),FEsolving_execIP(2)
crystallite_stress(ip,el) = all(crystallite_converged(:,ip,el))
enddo
enddo elementLooping5
@ -1129,27 +1129,27 @@ end function crystallite_stress
!--------------------------------------------------------------------------------------------------
!> @brief Backup data for homog cutback.
!--------------------------------------------------------------------------------------------------
subroutine constitutive_initializeRestorationPoints(i,e)
subroutine constitutive_initializeRestorationPoints(ip,el)
integer, intent(in) :: &
i, & !< integration point number
e !< element number
ip, & !< integration point number
el !< element number
integer :: &
c, & !< constituent number
co, & !< constituent number
s,ph, me
do c = 1,homogenization_Nconstituents(material_homogenizationAt(e))
ph = material_phaseAt(c,e)
me = material_phaseMemberAt(c,i,e)
crystallite_partitionedLp0(1:3,1:3,c,i,e) = crystallite_Lp0(1:3,1:3,c,i,e)
crystallite_partitionedF0(1:3,1:3,c,i,e) = crystallite_F0(1:3,1:3,c,i,e)
crystallite_partitionedS0(1:3,1:3,c,i,e) = crystallite_S0(1:3,1:3,c,i,e)
do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
crystallite_partitionedLp0(1:3,1:3,co,ip,el) = crystallite_Lp0(1:3,1:3,co,ip,el)
crystallite_partitionedF0(1:3,1:3,co,ip,el) = crystallite_F0(1:3,1:3,co,ip,el)
crystallite_partitionedS0(1:3,1:3,co,ip,el) = crystallite_S0(1:3,1:3,co,ip,el)
call mech_initializeRestorationPoints(ph,me)
do s = 1, phase_Nsources(material_phaseAt(c,e))
sourceState(material_phaseAt(c,e))%p(s)%partitionedState0(:,material_phasememberAt(c,i,e)) = &
sourceState(material_phaseAt(c,e))%p(s)%state0( :,material_phasememberAt(c,i,e))
do s = 1, phase_Nsources(material_phaseAt(co,el))
sourceState(material_phaseAt(co,el))%p(s)%partitionedState0(:,material_phasememberAt(co,ip,el)) = &
sourceState(material_phaseAt(co,el))%p(s)%state0( :,material_phasememberAt(co,ip,el))
enddo
enddo
@ -1159,23 +1159,23 @@ end subroutine constitutive_initializeRestorationPoints
!--------------------------------------------------------------------------------------------------
!> @brief Wind homog inc forward.
!--------------------------------------------------------------------------------------------------
subroutine constitutive_windForward(i,e)
subroutine constitutive_windForward(ip,el)
integer, intent(in) :: &
i, & !< integration point number
e !< element number
ip, & !< integration point number
el !< element number
integer :: &
c, & !< constituent number
co, & !< constituent number
s, ph, me
do c = 1,homogenization_Nconstituents(material_homogenizationAt(e))
ph = material_phaseAt(c,e)
me = material_phaseMemberAt(c,i,e)
crystallite_partitionedF0 (1:3,1:3,c,i,e) = crystallite_partitionedF(1:3,1:3,c,i,e)
crystallite_partitionedLp0(1:3,1:3,c,i,e) = crystallite_Lp (1:3,1:3,c,i,e)
crystallite_partitionedS0 (1:3,1:3,c,i,e) = crystallite_S (1:3,1:3,c,i,e)
do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
crystallite_partitionedF0 (1:3,1:3,co,ip,el) = crystallite_partitionedF(1:3,1:3,co,ip,el)
crystallite_partitionedLp0(1:3,1:3,co,ip,el) = crystallite_Lp (1:3,1:3,co,ip,el)
crystallite_partitionedS0 (1:3,1:3,co,ip,el) = crystallite_S (1:3,1:3,co,ip,el)
call constitutive_mech_windForward(ph,me)
do s = 1, phase_Nsources(material_phaseAt(c,e))
do s = 1, phase_Nsources(material_phaseAt(co,el))
sourceState(ph)%p(s)%partitionedState0(:,me) = sourceState(ph)%p(s)%state(:,me)
enddo
enddo
@ -1186,30 +1186,30 @@ end subroutine constitutive_windForward
!--------------------------------------------------------------------------------------------------
!> @brief Restore data after homog cutback.
!--------------------------------------------------------------------------------------------------
subroutine crystallite_restore(i,e,includeL)
subroutine crystallite_restore(ip,el,includeL)
integer, intent(in) :: &
i, & !< integration point number
e !< element number
ip, & !< integration point number
el !< element number
logical, intent(in) :: &
includeL !< protect agains fake cutback
integer :: &
c, p, m !< constituent number
co, p, m !< constituent number
do c = 1,homogenization_Nconstituents(material_homogenizationAt(e))
p = material_phaseAt(c,e)
m = material_phaseMemberAt(c,i,e)
do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
p = material_phaseAt(co,el)
m = material_phaseMemberAt(co,ip,el)
if (includeL) then
crystallite_Lp(1:3,1:3,c,i,e) = crystallite_partitionedLp0(1:3,1:3,c,i,e)
crystallite_Lp(1:3,1:3,co,ip,el) = crystallite_partitionedLp0(1:3,1:3,co,ip,el)
constitutive_mech_Li(p)%data(1:3,1:3,m) = constitutive_mech_partionedLi0(p)%data(1:3,1:3,m)
endif ! maybe protecting everything from overwriting makes more sense
constitutive_mech_Fp(p)%data(1:3,1:3,m) = constitutive_mech_partionedFp0(p)%data(1:3,1:3,m)
constitutive_mech_Fi(p)%data(1:3,1:3,m) = constitutive_mech_partionedFi0(p)%data(1:3,1:3,m)
crystallite_S (1:3,1:3,c,i,e) = crystallite_partitionedS0 (1:3,1:3,c,i,e)
crystallite_S (1:3,1:3,co,ip,el) = crystallite_partitionedS0 (1:3,1:3,co,ip,el)
plasticState (material_phaseAt(c,e))%state( :,material_phasememberAt(c,i,e)) = &
plasticState (material_phaseAt(c,e))%partitionedState0(:,material_phasememberAt(c,i,e))
plasticState (material_phaseAt(co,el))%state( :,material_phasememberAt(co,ip,el)) = &
plasticState (material_phaseAt(co,el))%partitionedState0(:,material_phasememberAt(co,ip,el))
enddo
end subroutine crystallite_restore
@ -1218,13 +1218,13 @@ end subroutine crystallite_restore
!--------------------------------------------------------------------------------------------------
!> @brief Calculate tangent (dPdF).
!--------------------------------------------------------------------------------------------------
function crystallite_stressTangent(c,i,e) result(dPdF)
function crystallite_stressTangent(co,ip,el) result(dPdF)
real(pReal), dimension(3,3,3,3) :: dPdF
integer, intent(in) :: &
c, & !< counter in constituent loop
i, & !< counter in integration point loop
e !< counter in element loop
co, & !< counter in constituent loop
ip, & !< counter in integration point loop
el !< counter in element loop
integer :: &
o, &
p, pp, m
@ -1247,21 +1247,21 @@ function crystallite_stressTangent(c,i,e) result(dPdF)
real(pReal), dimension(9,9):: temp_99
logical :: error
pp = material_phaseAt(c,e)
m = material_phaseMemberAt(c,i,e)
pp = material_phaseAt(co,el)
m = material_phaseMemberAt(co,ip,el)
call constitutive_hooke_SandItsTangents(devNull,dSdFe,dSdFi, &
crystallite_Fe(1:3,1:3,c,i,e), &
constitutive_mech_Fi(pp)%data(1:3,1:3,m),c,i,e)
crystallite_Fe(1:3,1:3,co,ip,el), &
constitutive_mech_Fi(pp)%data(1:3,1:3,m),co,ip,el)
call constitutive_LiAndItsTangents(devNull,dLidS,dLidFi, &
crystallite_S (1:3,1:3,c,i,e), &
crystallite_S (1:3,1:3,co,ip,el), &
constitutive_mech_Fi(pp)%data(1:3,1:3,m), &
c,i,e)
co,ip,el)
invFp = math_inv33(constitutive_mech_Fp(pp)%data(1:3,1:3,m))
invFi = math_inv33(constitutive_mech_Fi(pp)%data(1:3,1:3,m))
invSubFp0 = math_inv33(crystallite_subFp0(1:3,1:3,c,i,e))
invSubFi0 = math_inv33(crystallite_subFi0(1:3,1:3,c,i,e))
invSubFp0 = math_inv33(crystallite_subFp0(1:3,1:3,co,ip,el))
invSubFi0 = math_inv33(crystallite_subFi0(1:3,1:3,co,ip,el))
if (sum(abs(dLidS)) < tol_math_check) then
dFidS = 0.0_pReal
@ -1269,15 +1269,15 @@ function crystallite_stressTangent(c,i,e) result(dPdF)
lhs_3333 = 0.0_pReal; rhs_3333 = 0.0_pReal
do o=1,3; do p=1,3
lhs_3333(1:3,1:3,o,p) = lhs_3333(1:3,1:3,o,p) &
+ crystallite_subdt(c,i,e)*matmul(invSubFi0,dLidFi(1:3,1:3,o,p))
+ crystallite_subdt(co,ip,el)*matmul(invSubFi0,dLidFi(1:3,1:3,o,p))
lhs_3333(1:3,o,1:3,p) = lhs_3333(1:3,o,1:3,p) &
+ invFi*invFi(p,o)
rhs_3333(1:3,1:3,o,p) = rhs_3333(1:3,1:3,o,p) &
- crystallite_subdt(c,i,e)*matmul(invSubFi0,dLidS(1:3,1:3,o,p))
- crystallite_subdt(co,ip,el)*matmul(invSubFi0,dLidS(1:3,1:3,o,p))
enddo; enddo
call math_invert(temp_99,error,math_3333to99(lhs_3333))
if (error) then
call IO_warning(warning_ID=600,el=e,ip=i,g=c, &
call IO_warning(warning_ID=600,el=el,ip=ip,g=co, &
ext_msg='inversion error in analytic tangent calculation')
dFidS = 0.0_pReal
else
@ -1287,27 +1287,27 @@ function crystallite_stressTangent(c,i,e) result(dPdF)
endif
call constitutive_plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, &
crystallite_S (1:3,1:3,c,i,e), &
constitutive_mech_Fi(pp)%data(1:3,1:3,m),c,i,e)
crystallite_S (1:3,1:3,co,ip,el), &
constitutive_mech_Fi(pp)%data(1:3,1:3,m),co,ip,el)
dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS
!--------------------------------------------------------------------------------------------------
! calculate dSdF
temp_33_1 = transpose(matmul(invFp,invFi))
temp_33_2 = matmul(crystallite_subF(1:3,1:3,c,i,e),invSubFp0)
temp_33_3 = matmul(matmul(crystallite_subF(1:3,1:3,c,i,e),invFp), invSubFi0)
temp_33_2 = matmul(crystallite_subF(1:3,1:3,co,ip,el),invSubFp0)
temp_33_3 = matmul(matmul(crystallite_subF(1:3,1:3,co,ip,el),invFp), invSubFi0)
do o=1,3; do p=1,3
rhs_3333(p,o,1:3,1:3) = matmul(dSdFe(p,o,1:3,1:3),temp_33_1)
temp_3333(1:3,1:3,p,o) = matmul(matmul(temp_33_2,dLpdS(1:3,1:3,p,o)), invFi) &
+ matmul(temp_33_3,dLidS(1:3,1:3,p,o))
enddo; enddo
lhs_3333 = crystallite_subdt(c,i,e)*math_mul3333xx3333(dSdFe,temp_3333) &
lhs_3333 = crystallite_subdt(co,ip,el)*math_mul3333xx3333(dSdFe,temp_3333) &
+ math_mul3333xx3333(dSdFi,dFidS)
call math_invert(temp_99,error,math_eye(9)+math_3333to99(lhs_3333))
if (error) then
call IO_warning(warning_ID=600,el=e,ip=i,g=c, &
call IO_warning(warning_ID=600,el=el,ip=ip,g=co, &
ext_msg='inversion error in analytic tangent calculation')
dSdF = rhs_3333
else
@ -1318,16 +1318,16 @@ function crystallite_stressTangent(c,i,e) result(dPdF)
! calculate dFpinvdF
temp_3333 = math_mul3333xx3333(dLpdS,dSdF)
do o=1,3; do p=1,3
dFpinvdF(1:3,1:3,p,o) = -crystallite_subdt(c,i,e) &
dFpinvdF(1:3,1:3,p,o) = -crystallite_subdt(co,ip,el) &
* matmul(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi))
enddo; enddo
!--------------------------------------------------------------------------------------------------
! assemble dPdF
temp_33_1 = matmul(crystallite_S(1:3,1:3,c,i,e),transpose(invFp))
temp_33_1 = matmul(crystallite_S(1:3,1:3,co,ip,el),transpose(invFp))
temp_33_2 = matmul(invFp,temp_33_1)
temp_33_3 = matmul(crystallite_subF(1:3,1:3,c,i,e),invFp)
temp_33_4 = matmul(temp_33_3,crystallite_S(1:3,1:3,c,i,e))
temp_33_3 = matmul(crystallite_subF(1:3,1:3,co,ip,el),invFp)
temp_33_4 = matmul(temp_33_3,crystallite_S(1:3,1:3,co,ip,el))
dPdF = 0.0_pReal
do p=1,3
@ -1335,7 +1335,7 @@ function crystallite_stressTangent(c,i,e) result(dPdF)
enddo
do o=1,3; do p=1,3
dPdF(1:3,1:3,p,o) = dPdF(1:3,1:3,p,o) &
+ matmul(matmul(crystallite_subF(1:3,1:3,c,i,e), &
+ matmul(matmul(crystallite_subF(1:3,1:3,co,ip,el), &
dFpinvdF(1:3,1:3,p,o)),temp_33_1) &
+ matmul(matmul(temp_33_3,dSdF(1:3,1:3,p,o)), &
transpose(invFp)) &
@ -1351,25 +1351,26 @@ end function crystallite_stressTangent
subroutine crystallite_orientations
integer &
c, & !< counter in integration point component loop
i, & !< counter in integration point loop
e !< counter in element loop
co, & !< counter in integration point component loop
ip, & !< counter in integration point loop
el !< counter in element loop
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do c = 1,homogenization_Nconstituents(material_homogenizationAt(e))
call crystallite_orientation(c,i,e)%fromMatrix(transpose(math_rotationalPart(crystallite_Fe(1:3,1:3,c,i,e))))
do el = FEsolving_execElem(1),FEsolving_execElem(2)
do ip = FEsolving_execIP(1),FEsolving_execIP(2)
do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
call crystallite_orientation(co,ip,el)%fromMatrix(transpose(math_rotationalPart(crystallite_Fe(1:3,1:3,co,ip,el))))
enddo; enddo; enddo
!$OMP END PARALLEL DO
nonlocalPresent: if (any(plasticState%nonlocal)) then
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2)
if (plasticState(material_phaseAt(1,e))%nonlocal) then
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do el = FEsolving_execElem(1),FEsolving_execElem(2)
if (plasticState(material_phaseAt(1,el))%nonlocal) then
do ip = FEsolving_execIP(1),FEsolving_execIP(2)
call plastic_nonlocal_updateCompatibility(crystallite_orientation, &
phase_plasticityInstance(material_phaseAt(1,e)),i,e)
phase_plasticityInstance(material_phaseAt(1,el)),ip,el)
enddo
endif
enddo
@ -1403,15 +1404,15 @@ end function crystallite_push33ToRef
!> @brief integrate stress, state with adaptive 1st order explicit Euler method
!> using Fixed Point Iteration to adapt the stepsize
!--------------------------------------------------------------------------------------------------
subroutine integrateSourceState(g,i,e)
subroutine integrateSourceState(co,ip,el)
integer, intent(in) :: &
e, & !< element index in element loop
i, & !< integration point index in ip loop
g !< grain index in grain loop
el, & !< element index in element loop
ip, & !< integration point index in ip loop
co !< grain index in grain loop
integer :: &
NiterationState, & !< number of iterations in state loop
p, &
ph, &
c, &
s, &
size_pl
@ -1425,51 +1426,51 @@ subroutine integrateSourceState(g,i,e)
logical :: &
broken
p = material_phaseAt(g,e)
c = material_phaseMemberAt(g,i,e)
ph = material_phaseAt(co,el)
c = material_phaseMemberAt(co,ip,el)
broken = constitutive_thermal_collectDotState(p,c)
broken = broken .or. constitutive_damage_collectDotState(crystallite_S(1:3,1:3,g,i,e), g,i,e,p,c)
broken = constitutive_thermal_collectDotState(ph,c)
broken = broken .or. constitutive_damage_collectDotState(crystallite_S(1:3,1:3,co,ip,el), co,ip,el,ph,c)
if(broken) return
do s = 1, phase_Nsources(p)
size_so(s) = sourceState(p)%p(s)%sizeDotState
sourceState(p)%p(s)%state(1:size_so(s),c) = sourceState(p)%p(s)%subState0(1:size_so(s),c) &
+ sourceState(p)%p(s)%dotState (1:size_so(s),c) &
* crystallite_subdt(g,i,e)
do s = 1, phase_Nsources(ph)
size_so(s) = sourceState(ph)%p(s)%sizeDotState
sourceState(ph)%p(s)%state(1:size_so(s),c) = sourceState(ph)%p(s)%subState0(1:size_so(s),c) &
+ sourceState(ph)%p(s)%dotState (1:size_so(s),c) &
* crystallite_subdt(co,ip,el)
source_dotState(1:size_so(s),2,s) = 0.0_pReal
enddo
iteration: do NiterationState = 1, num%nState
do s = 1, phase_Nsources(p)
do s = 1, phase_Nsources(ph)
if(nIterationState > 1) source_dotState(1:size_so(s),2,s) = source_dotState(1:size_so(s),1,s)
source_dotState(1:size_so(s),1,s) = sourceState(p)%p(s)%dotState(:,c)
source_dotState(1:size_so(s),1,s) = sourceState(ph)%p(s)%dotState(:,c)
enddo
broken = constitutive_thermal_collectDotState(p,c)
broken = broken .or. constitutive_damage_collectDotState(crystallite_S(1:3,1:3,g,i,e), g,i,e,p,c)
broken = constitutive_thermal_collectDotState(ph,c)
broken = broken .or. constitutive_damage_collectDotState(crystallite_S(1:3,1:3,co,ip,el), co,ip,el,ph,c)
if(broken) exit iteration
do s = 1, phase_Nsources(p)
zeta = damper(sourceState(p)%p(s)%dotState(:,c), &
do s = 1, phase_Nsources(ph)
zeta = damper(sourceState(ph)%p(s)%dotState(:,c), &
source_dotState(1:size_so(s),1,s),&
source_dotState(1:size_so(s),2,s))
sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) * zeta &
sourceState(ph)%p(s)%dotState(:,c) = sourceState(ph)%p(s)%dotState(:,c) * zeta &
+ source_dotState(1:size_so(s),1,s)* (1.0_pReal - zeta)
r(1:size_so(s)) = sourceState(p)%p(s)%state (1:size_so(s),c) &
- sourceState(p)%p(s)%subState0(1:size_so(s),c) &
- sourceState(p)%p(s)%dotState (1:size_so(s),c) * crystallite_subdt(g,i,e)
sourceState(p)%p(s)%state(1:size_so(s),c) = sourceState(p)%p(s)%state(1:size_so(s),c) &
r(1:size_so(s)) = sourceState(ph)%p(s)%state (1:size_so(s),c) &
- sourceState(ph)%p(s)%subState0(1:size_so(s),c) &
- sourceState(ph)%p(s)%dotState (1:size_so(s),c) * crystallite_subdt(co,ip,el)
sourceState(ph)%p(s)%state(1:size_so(s),c) = sourceState(ph)%p(s)%state(1:size_so(s),c) &
- r(1:size_so(s))
crystallite_converged(g,i,e) = &
crystallite_converged(g,i,e) .and. converged(r(1:size_so(s)), &
sourceState(p)%p(s)%state(1:size_so(s),c), &
sourceState(p)%p(s)%atol(1:size_so(s)))
crystallite_converged(co,ip,el) = &
crystallite_converged(co,ip,el) .and. converged(r(1:size_so(s)), &
sourceState(ph)%p(s)%state(1:size_so(s),c), &
sourceState(ph)%p(s)%atol(1:size_so(s)))
enddo
if(crystallite_converged(g,i,e)) then
broken = constitutive_damage_deltaState(crystallite_Fe(1:3,1:3,g,i,e),g,i,e,p,c)
if(crystallite_converged(co,ip,el)) then
broken = constitutive_damage_deltaState(crystallite_Fe(1:3,1:3,co,ip,el),co,ip,el,ph,c)
exit iteration
endif

288
src/constitutive_mech.f90
View File

@ -284,7 +284,7 @@ contains
module subroutine mech_init
integer :: &
p, &
ph, &
stiffDegradationCtr
class(tNode), pointer :: &
num_crystallite, &
@ -304,35 +304,35 @@ module subroutine mech_init
allocate(phase_NstiffnessDegradations(phases%length),source=0)
allocate(output_constituent(phases%length))
do p = 1, phases%length
phase => phases%get(p)
do ph = 1, phases%length
phase => phases%get(ph)
mech => phase%get('mechanics')
#if defined(__GFORTRAN__)
output_constituent(p)%label = output_asStrings(mech)
output_constituent(ph)%label = output_asStrings(mech)
#else
output_constituent(p)%label = mech%get_asStrings('output',defaultVal=emptyStringArray)
output_constituent(ph)%label = mech%get_asStrings('output',defaultVal=emptyStringArray)
#endif
elastic => mech%get('elasticity')
if(elastic%get_asString('type') == 'hooke') then
phase_elasticity(p) = ELASTICITY_HOOKE_ID
phase_elasticity(ph) = ELASTICITY_HOOKE_ID
else
call IO_error(200,ext_msg=elastic%get_asString('type'))
endif
stiffDegradation => mech%get('stiffness_degradation',defaultVal=emptyList) ! check for stiffness degradation mechanisms
phase_NstiffnessDegradations(p) = stiffDegradation%length
phase_NstiffnessDegradations(ph) = stiffDegradation%length
enddo
allocate(phase_stiffnessDegradation(maxval(phase_NstiffnessDegradations),phases%length), &
source=STIFFNESS_DEGRADATION_undefined_ID)
if(maxVal(phase_NstiffnessDegradations)/=0) then
do p = 1, phases%length
phase => phases%get(p)
do ph = 1, phases%length
phase => phases%get(ph)
mech => phase%get('mechanics')
stiffDegradation => mech%get('stiffness_degradation',defaultVal=emptyList)
do stiffDegradationCtr = 1, stiffDegradation%length
if(stiffDegradation%get_asString(stiffDegradationCtr) == 'damage') &
phase_stiffnessDegradation(stiffDegradationCtr,p) = STIFFNESS_DEGRADATION_damage_ID
phase_stiffnessDegradation(stiffDegradationCtr,ph) = STIFFNESS_DEGRADATION_damage_ID
enddo
enddo
endif
@ -352,9 +352,9 @@ module subroutine mech_init
where(plastic_dislotungsten_init()) phase_plasticity = PLASTICITY_DISLOTUNGSTEN_ID
where(plastic_nonlocal_init()) phase_plasticity = PLASTICITY_NONLOCAL_ID
do p = 1, phases%length
phase_elasticityInstance(p) = count(phase_elasticity(1:p) == phase_elasticity(p))
phase_plasticityInstance(p) = count(phase_plasticity(1:p) == phase_plasticity(p))
do ph = 1, phases%length
phase_elasticityInstance(ph) = count(phase_elasticity(1:ph) == phase_elasticity(ph))
phase_plasticityInstance(ph) = count(phase_plasticity(1:ph) == phase_plasticity(ph))
enddo
num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict)
@ -397,15 +397,15 @@ function plastic_active(plastic_label) result(active_plastic)
phase, &
mech, &
pl
integer :: p
integer :: ph
phases => config_material%get('phase')
allocate(active_plastic(phases%length), source = .false. )
do p = 1, phases%length
phase => phases%get(p)
do ph = 1, phases%length
phase => phases%get(ph)
mech => phase%get('mechanics')
pl => mech%get('plasticity')
if(pl%get_asString('type') == plastic_label) active_plastic(p) = .true.
if(pl%get_asString('type') == plastic_label) active_plastic(ph) = .true.
enddo
end function plastic_active
@ -568,13 +568,13 @@ end subroutine constitutive_plastic_LpAndItsTangents
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
function mech_collectDotState(subdt, co, ip, el,phase,of) result(broken)
function mech_collectDotState(subdt, co, ip, el,ph,of) result(broken)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el, & !< element
phase, &
ph, &
of
real(pReal), intent(in) :: &
subdt !< timestep
@ -588,12 +588,12 @@ function mech_collectDotState(subdt, co, ip, el,phase,of) result(broken)
logical :: broken
ho = material_homogenizationAt(el)
tme = material_homogenizationMemberAt(ip,el)
instance = phase_plasticityInstance(phase)
instance = phase_plasticityInstance(ph)
Mp = matmul(matmul(transpose(constitutive_mech_Fi(phase)%data(1:3,1:3,of)),&
constitutive_mech_Fi(phase)%data(1:3,1:3,of)),crystallite_S(1:3,1:3,co,ip,el))
Mp = matmul(matmul(transpose(constitutive_mech_Fi(ph)%data(1:3,1:3,of)),&
constitutive_mech_Fi(ph)%data(1:3,1:3,of)),crystallite_S(1:3,1:3,co,ip,el))
plasticityType: select case (phase_plasticity(phase))
plasticityType: select case (phase_plasticity(ph))
case (PLASTICITY_ISOTROPIC_ID) plasticityType
call plastic_isotropic_dotState(Mp,instance,of)
@ -614,7 +614,7 @@ function mech_collectDotState(subdt, co, ip, el,phase,of) result(broken)
call plastic_nonlocal_dotState(Mp,crystallite_partitionedF0,temperature(ho)%p(tme),subdt, &
instance,of,ip,el)
end select plasticityType
broken = any(IEEE_is_NaN(plasticState(phase)%dotState(:,of)))
broken = any(IEEE_is_NaN(plasticState(ph)%dotState(:,of)))
end function mech_collectDotState
@ -624,13 +624,13 @@ end function mech_collectDotState
!> @brief for constitutive models having an instantaneous change of state
!> will return false if delta state is not needed/supported by the constitutive model
!--------------------------------------------------------------------------------------------------
function constitutive_deltaState(S, Fi, co, ip, el, phase, of) result(broken)
function constitutive_deltaState(S, Fi, co, ip, el, ph, of) result(broken)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el, & !< element
phase, &
ph, &
of
real(pReal), intent(in), dimension(3,3) :: &
S, & !< 2nd Piola Kirchhoff stress
@ -645,17 +645,17 @@ function constitutive_deltaState(S, Fi, co, ip, el, phase, of) result(broken)
broken
Mp = matmul(matmul(transpose(Fi),Fi),S)
instance = phase_plasticityInstance(phase)
instance = phase_plasticityInstance(ph)
plasticityType: select case (phase_plasticity(phase))
plasticityType: select case (phase_plasticity(ph))
case (PLASTICITY_KINEHARDENING_ID) plasticityType
call plastic_kinehardening_deltaState(Mp,instance,of)
broken = any(IEEE_is_NaN(plasticState(phase)%deltaState(:,of)))
broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,of)))
case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_deltaState(Mp,instance,of,ip,el)
broken = any(IEEE_is_NaN(plasticState(phase)%deltaState(:,of)))
broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,of)))
case default
broken = .false.
@ -663,13 +663,13 @@ function constitutive_deltaState(S, Fi, co, ip, el, phase, of) result(broken)
end select plasticityType
if(.not. broken) then
select case(phase_plasticity(phase))
select case(phase_plasticity(ph))
case (PLASTICITY_NONLOCAL_ID,PLASTICITY_KINEHARDENING_ID)
myOffset = plasticState(phase)%offsetDeltaState
mySize = plasticState(phase)%sizeDeltaState
plasticState(phase)%state(myOffset + 1:myOffset + mySize,of) = &
plasticState(phase)%state(myOffset + 1:myOffset + mySize,of) + plasticState(phase)%deltaState(1:mySize,of)
myOffset = plasticState(ph)%offsetDeltaState
mySize = plasticState(ph)%sizeDeltaState
plasticState(ph)%state(myOffset + 1:myOffset + mySize,of) = &
plasticState(ph)%state(myOffset + 1:myOffset + mySize,of) + plasticState(ph)%deltaState(1:mySize,of)
end select
endif
@ -944,16 +944,16 @@ end function integrateStress
!> @brief integrate stress, state with adaptive 1st order explicit Euler method
!> using Fixed Point Iteration to adapt the stepsize
!--------------------------------------------------------------------------------------------------
module subroutine integrateStateFPI(g,i,e)
module subroutine integrateStateFPI(co,ip,el)
integer, intent(in) :: &
e, & !< element index in element loop
i, & !< integration point index in ip loop
g !< grain index in grain loop
el, & !< element index in element loop
ip, & !< integration point index in ip loop
co !< grain index in grain loop
integer :: &
NiterationState, & !< number of iterations in state loop
p, &
c, &
ph, &
me, &
s, &
size_pl
integer, dimension(maxval(phase_Nsources)) :: &
@ -968,45 +968,45 @@ module subroutine integrateStateFPI(g,i,e)
logical :: &
broken
p = material_phaseAt(g,e)
c = material_phaseMemberAt(g,i,e)
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
broken = mech_collectDotState(crystallite_subdt(g,i,e), g,i,e,p,c)
broken = mech_collectDotState(crystallite_subdt(co,ip,el), co,ip,el,ph,me)
if(broken) return
size_pl = plasticState(p)%sizeDotState
plasticState(p)%state(1:size_pl,c) = plasticState(p)%subState0(1:size_pl,c) &
+ plasticState(p)%dotState (1:size_pl,c) &
* crystallite_subdt(g,i,e)
size_pl = plasticState(ph)%sizeDotState
plasticState(ph)%state(1:size_pl,me) = plasticState(ph)%subState0(1:size_pl,me) &
+ plasticState(ph)%dotState (1:size_pl,me) &
* crystallite_subdt(co,ip,el)
plastic_dotState(1:size_pl,2) = 0.0_pReal
iteration: do NiterationState = 1, num%nState
if(nIterationState > 1) plastic_dotState(1:size_pl,2) = plastic_dotState(1:size_pl,1)
plastic_dotState(1:size_pl,1) = plasticState(p)%dotState(:,c)
plastic_dotState(1:size_pl,1) = plasticState(ph)%dotState(:,me)
broken = integrateStress(g,i,e)
broken = integrateStress(co,ip,el)
if(broken) exit iteration
broken = mech_collectDotState(crystallite_subdt(g,i,e), g,i,e,p,c)
broken = mech_collectDotState(crystallite_subdt(co,ip,el), co,ip,el,ph,me)
if(broken) exit iteration
zeta = damper(plasticState(p)%dotState(:,c),plastic_dotState(1:size_pl,1),&
zeta = damper(plasticState(ph)%dotState(:,me),plastic_dotState(1:size_pl,1),&
plastic_dotState(1:size_pl,2))
plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * zeta &
plasticState(ph)%dotState(:,me) = plasticState(ph)%dotState(:,me) * zeta &
+ plastic_dotState(1:size_pl,1) * (1.0_pReal - zeta)
r(1:size_pl) = plasticState(p)%state (1:size_pl,c) &
- plasticState(p)%subState0(1:size_pl,c) &
- plasticState(p)%dotState (1:size_pl,c) * crystallite_subdt(g,i,e)
plasticState(p)%state(1:size_pl,c) = plasticState(p)%state(1:size_pl,c) &
r(1:size_pl) = plasticState(ph)%state (1:size_pl,me) &
- plasticState(ph)%subState0(1:size_pl,me) &
- plasticState(ph)%dotState (1:size_pl,me) * crystallite_subdt(co,ip,el)
plasticState(ph)%state(1:size_pl,me) = plasticState(ph)%state(1:size_pl,me) &
- r(1:size_pl)
crystallite_converged(g,i,e) = converged(r(1:size_pl), &
plasticState(p)%state(1:size_pl,c), &
plasticState(p)%atol(1:size_pl))
crystallite_converged(co,ip,el) = converged(r(1:size_pl), &
plasticState(ph)%state(1:size_pl,me), &
plasticState(ph)%atol(1:size_pl))
if(crystallite_converged(g,i,e)) then
broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), &
constitutive_mech_Fi(p)%data(1:3,1:3,c),g,i,e,p,c)
if(crystallite_converged(co,ip,el)) then
broken = constitutive_deltaState(crystallite_S(1:3,1:3,co,ip,el), &
constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el,ph,me)
exit iteration
endif
@ -1041,36 +1041,36 @@ end subroutine integrateStateFPI
!--------------------------------------------------------------------------------------------------
!> @brief integrate state with 1st order explicit Euler method
!--------------------------------------------------------------------------------------------------
subroutine integrateStateEuler(g,i,e)
subroutine integrateStateEuler(co,ip,el)
integer, intent(in) :: &
e, & !< element index in element loop
i, & !< integration point index in ip loop
g !< grain index in grain loop
el, & !< element index in element loop
ip, & !< integration point index in ip loop
co !< grain index in grain loop
integer :: &
p, &
c, &
ph, &
me, &
sizeDotState
logical :: &
broken
p = material_phaseAt(g,e)
c = material_phaseMemberAt(g,i,e)
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
broken = mech_collectDotState(crystallite_subdt(g,i,e), g,i,e,p,c)
broken = mech_collectDotState(crystallite_subdt(co,ip,el), co,ip,el,ph,me)
if(broken) return
sizeDotState = plasticState(p)%sizeDotState
plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) &
+ plasticState(p)%dotState (1:sizeDotState,c) &
* crystallite_subdt(g,i,e)
sizeDotState = plasticState(ph)%sizeDotState
plasticState(ph)%state(1:sizeDotState,me) = plasticState(ph)%subState0(1:sizeDotState,me) &
+ plasticState(ph)%dotState (1:sizeDotState,me) &
* crystallite_subdt(co,ip,el)
broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), &
constitutive_mech_Fi(p)%data(1:3,1:3,c),g,i,e,p,c)
broken = constitutive_deltaState(crystallite_S(1:3,1:3,co,ip,el), &
constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el,ph,me)
if(broken) return
broken = integrateStress(g,i,e)
crystallite_converged(g,i,e) = .not. broken
broken = integrateStress(co,ip,el)
crystallite_converged(co,ip,el) = .not. broken
end subroutine integrateStateEuler
@ -1078,15 +1078,15 @@ end subroutine integrateStateEuler
!--------------------------------------------------------------------------------------------------
!> @brief integrate stress, state with 1st order Euler method with adaptive step size
!--------------------------------------------------------------------------------------------------
subroutine integrateStateAdaptiveEuler(g,i,e)
subroutine integrateStateAdaptiveEuler(co,ip,el)
integer, intent(in) :: &
e, & !< element index in element loop
i, & !< integration point index in ip loop
g !< grain index in grain loop
el, & !< element index in element loop
ip, & !< integration point index in ip loop
co !< grain index in grain loop
integer :: &
p, &
c, &
ph, &
me, &
sizeDotState
logical :: &
broken
@ -1094,34 +1094,34 @@ subroutine integrateStateAdaptiveEuler(g,i,e)
real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: residuum_plastic
p = material_phaseAt(g,e)
c = material_phaseMemberAt(g,i,e)
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
broken = mech_collectDotState(crystallite_subdt(g,i,e), g,i,e,p,c)
broken = mech_collectDotState(crystallite_subdt(co,ip,el), co,ip,el,ph,me)
if(broken) return
sizeDotState = plasticState(p)%sizeDotState
sizeDotState = plasticState(ph)%sizeDotState
residuum_plastic(1:sizeDotState) = - plasticState(p)%dotstate(1:sizeDotState,c) * 0.5_pReal * crystallite_subdt(g,i,e)
plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) &
+ plasticState(p)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e)
residuum_plastic(1:sizeDotState) = - plasticState(ph)%dotstate(1:sizeDotState,me) * 0.5_pReal * crystallite_subdt(co,ip,el)
plasticState(ph)%state(1:sizeDotState,me) = plasticState(ph)%subState0(1:sizeDotState,me) &
+ plasticState(ph)%dotstate(1:sizeDotState,me) * crystallite_subdt(co,ip,el)
broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), &
constitutive_mech_Fi(p)%data(1:3,1:3,c),g,i,e,p,c)
broken = constitutive_deltaState(crystallite_S(1:3,1:3,co,ip,el), &
constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el,ph,me)
if(broken) return
broken = integrateStress(g,i,e)
broken = integrateStress(co,ip,el)
if(broken) return
broken = mech_collectDotState(crystallite_subdt(g,i,e), g,i,e,p,c)
broken = mech_collectDotState(crystallite_subdt(co,ip,el), co,ip,el,ph,me)
if(broken) return
sizeDotState = plasticState(p)%sizeDotState
crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState) &
+ 0.5_pReal * plasticState(p)%dotState(:,c) * crystallite_subdt(g,i,e), &
plasticState(p)%state(1:sizeDotState,c), &
plasticState(p)%atol(1:sizeDotState))
sizeDotState = plasticState(ph)%sizeDotState
crystallite_converged(co,ip,el) = converged(residuum_plastic(1:sizeDotState) &
+ 0.5_pReal * plasticState(ph)%dotState(:,me) * crystallite_subdt(co,ip,el), &
plasticState(ph)%state(1:sizeDotState,me), &
plasticState(ph)%atol(1:sizeDotState))
end subroutine integrateStateAdaptiveEuler
@ -1129,9 +1129,9 @@ end subroutine integrateStateAdaptiveEuler
!---------------------------------------------------------------------------------------------------
!> @brief Integrate state (including stress integration) with the classic Runge Kutta method
!---------------------------------------------------------------------------------------------------
subroutine integrateStateRK4(g,i,e)
subroutine integrateStateRK4(co,ip,el)
integer, intent(in) :: g,i,e
integer, intent(in) :: co,ip,el
real(pReal), dimension(3,3), parameter :: &
A = reshape([&
@ -1144,7 +1144,7 @@ subroutine integrateStateRK4(g,i,e)
real(pReal), dimension(4), parameter :: &
B = [1.0_pReal/6.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/6.0_pReal]
call integrateStateRK(g,i,e,A,B,C)
call integrateStateRK(co,ip,el,A,B,C)
end subroutine integrateStateRK4
@ -1152,9 +1152,9 @@ end subroutine integrateStateRK4
!---------------------------------------------------------------------------------------------------
!> @brief Integrate state (including stress integration) with the Cash-Carp method
!---------------------------------------------------------------------------------------------------
subroutine integrateStateRKCK45(g,i,e)
subroutine integrateStateRKCK45(co,ip,el)
integer, intent(in) :: g,i,e
integer, intent(in) :: co,ip,el
real(pReal), dimension(5,5), parameter :: &
A = reshape([&
@ -1174,7 +1174,7 @@ subroutine integrateStateRKCK45(g,i,e)
[2825.0_pReal/27648.0_pReal, .0_pReal, 18575.0_pReal/48384.0_pReal,&
13525.0_pReal/55296.0_pReal, 277.0_pReal/14336.0_pReal, 1._pReal/4._pReal]
call integrateStateRK(g,i,e,A,B,C,DB)
call integrateStateRK(co,ip,el,A,B,C,DB)
end subroutine integrateStateRKCK45
@ -1183,7 +1183,7 @@ end subroutine integrateStateRKCK45
!> @brief Integrate state (including stress integration) with an explicit Runge-Kutta method or an
!! embedded explicit Runge-Kutta method
!--------------------------------------------------------------------------------------------------
subroutine integrateStateRK(g,i,e,A,B,CC,DB)
subroutine integrateStateRK(co,ip,el,A,B,CC,DB)
real(pReal), dimension(:,:), intent(in) :: A
@ -1191,71 +1191,71 @@ subroutine integrateStateRK(g,i,e,A,B,CC,DB)
real(pReal), dimension(:), intent(in), optional :: DB
integer, intent(in) :: &
e, & !< element index in element loop
i, & !< integration point index in ip loop
g !< grain index in grain loop
el, & !< element index in element loop
ip, & !< integration point index in ip loop
co !< grain index in grain loop
integer :: &
stage, & ! stage index in integration stage loop
n, &
p, &
c, &
ph, &
me, &
sizeDotState
logical :: &
broken
real(pReal), dimension(constitutive_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState
p = material_phaseAt(g,e)
c = material_phaseMemberAt(g,i,e)
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
broken = mech_collectDotState(crystallite_subdt(g,i,e), g,i,e,p,c)
broken = mech_collectDotState(crystallite_subdt(co,ip,el), co,ip,el,ph,me)
if(broken) return
do stage = 1,size(A,1)
sizeDotState = plasticState(p)%sizeDotState
plastic_RKdotState(1:sizeDotState,stage) = plasticState(p)%dotState(:,c)
plasticState(p)%dotState(:,c) = A(1,stage) * plastic_RKdotState(1:sizeDotState,1)
sizeDotState = plasticState(ph)%sizeDotState
plastic_RKdotState(1:sizeDotState,stage) = plasticState(ph)%dotState(:,me)
plasticState(ph)%dotState(:,me) = A(1,stage) * plastic_RKdotState(1:sizeDotState,1)
do n = 2, stage
sizeDotState = plasticState(p)%sizeDotState
plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) &
sizeDotState = plasticState(ph)%sizeDotState
plasticState(ph)%dotState(:,me) = plasticState(ph)%dotState(:,me) &
+ A(n,stage) * plastic_RKdotState(1:sizeDotState,n)
enddo
sizeDotState = plasticState(p)%sizeDotState
plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) &
+ plasticState(p)%dotState (1:sizeDotState,c) &
* crystallite_subdt(g,i,e)
sizeDotState = plasticState(ph)%sizeDotState
plasticState(ph)%state(1:sizeDotState,me) = plasticState(ph)%subState0(1:sizeDotState,me) &
+ plasticState(ph)%dotState (1:sizeDotState,me) &
* crystallite_subdt(co,ip,el)
broken = integrateStress(g,i,e,CC(stage))
broken = integrateStress(co,ip,el,CC(stage))
if(broken) exit
broken = mech_collectDotState(crystallite_subdt(g,i,e)*CC(stage), g,i,e,p,c)
broken = mech_collectDotState(crystallite_subdt(co,ip,el)*CC(stage), co,ip,el,ph,me)
if(broken) exit
enddo
if(broken) return
sizeDotState = plasticState(p)%sizeDotState
sizeDotState = plasticState(ph)%sizeDotState
plastic_RKdotState(1:sizeDotState,size(B)) = plasticState (p)%dotState(:,c)
plasticState(p)%dotState(:,c) = matmul(plastic_RKdotState(1:sizeDotState,1:size(B)),B)
plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) &
+ plasticState(p)%dotState (1:sizeDotState,c) &
* crystallite_subdt(g,i,e)
plastic_RKdotState(1:sizeDotState,size(B)) = plasticState (ph)%dotState(:,me)
plasticState(ph)%dotState(:,me) = matmul(plastic_RKdotState(1:sizeDotState,1:size(B)),B)
plasticState(ph)%state(1:sizeDotState,me) = plasticState(ph)%subState0(1:sizeDotState,me) &
+ plasticState(ph)%dotState (1:sizeDotState,me) &
* crystallite_subdt(co,ip,el)
if(present(DB)) &
broken = .not. converged( matmul(plastic_RKdotState(1:sizeDotState,1:size(DB)),DB) &
* crystallite_subdt(g,i,e), &
plasticState(p)%state(1:sizeDotState,c), &
plasticState(p)%atol(1:sizeDotState))
* crystallite_subdt(co,ip,el), &
plasticState(ph)%state(1:sizeDotState,me), &
plasticState(ph)%atol(1:sizeDotState))
if(broken) return
broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), &
constitutive_mech_Fi(p)%data(1:3,1:3,c),g,i,e,p,c)
broken = constitutive_deltaState(crystallite_S(1:3,1:3,co,ip,el), &
constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el,ph,me)
if(broken) return
broken = integrateStress(g,i,e)
crystallite_converged(g,i,e) = .not. broken
broken = integrateStress(co,ip,el)
crystallite_converged(co,ip,el) = .not. broken
end subroutine integrateStateRK
@ -1396,9 +1396,7 @@ end subroutine crystallite_results
!--------------------------------------------------------------------------------------------------
module subroutine mech_initializeRestorationPoints(ph,me)
integer, intent(in) :: &
ph, &
me
integer, intent(in) :: ph, me
constitutive_mech_partionedFi0(ph)%data(1:3,1:3,me) = constitutive_mech_Fi0(ph)%data(1:3,1:3,me)
@ -1416,13 +1414,13 @@ module subroutine constitutive_mech_windForward(ph,me)
integer, intent(in) :: ph, me
constitutive_mech_partionedFp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp(ph)%data(1:3,1:3,me)
constitutive_mech_partionedFi0(ph)%data(1:3,1:3,me) = constitutive_mech_Fi(ph)%data(1:3,1:3,me)
constitutive_mech_partionedLi0(ph)%data(1:3,1:3,me) = constitutive_mech_Li(ph)%data(1:3,1:3,me)
plasticState(ph)%partitionedState0(:,me) = plasticState(ph)%state(:,me)
end subroutine constitutive_mech_windForward