diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 4e29f3f96..461a67c9a 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -117,7 +117,7 @@ module constitutive type(tDebugOptions) :: debugCrystallite - procedure(integrateStateFPI), pointer :: integrateState + integer(kind(PLASTICITY_undefined_ID)), dimension(:), allocatable, public :: & phase_plasticity !< plasticity of each phase @@ -187,6 +187,12 @@ module constitutive ! == cleaned:end =================================================================================== + module function crystallite_stress(dt,co,ip,el) + real(pReal), intent(in) :: dt + integer, intent(in) :: co, ip, el + logical :: crystallite_stress + end function crystallite_stress + module function constitutive_homogenizedC(co,ip,el) result(C) integer, intent(in) :: co, ip, el real(pReal), dimension(6,6) :: C @@ -362,10 +368,6 @@ module constitutive dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient end subroutine constitutive_hooke_SandItsTangents - module subroutine integrateStateFPI(co,ip,el) - integer, intent(in) :: co, ip, el - end subroutine integrateStateFPI - end interface @@ -392,6 +394,7 @@ module constitutive crystallite_orientations, & crystallite_push33ToRef, & crystallite_restartWrite, & + integrateSourceState, & crystallite_restartRead, & constitutive_initializeRestorationPoints, & constitutive_windForward, & @@ -983,120 +986,7 @@ subroutine crystallite_init end subroutine crystallite_init -!-------------------------------------------------------------------------------------------------- -!> @brief calculate stress (P) -!-------------------------------------------------------------------------------------------------- -function crystallite_stress(dt,co,ip,el) - real(pReal), intent(in) :: dt - integer, intent(in) :: & - co, & - ip, & - el - - logical :: crystallite_stress - - real(pReal) :: & - formerSubStep - integer :: & - NiterationCrystallite, & ! number of iterations in crystallite loop - s, ph, me - logical :: todo - real(pReal) :: subFrac,subStep - real(pReal), dimension(3,3) :: & - subLp0, & !< plastic velocity grad at start of crystallite inc - subLi0 !< intermediate velocity grad at start of crystallite inc - - - ph = material_phaseAt(co,el) - me = material_phaseMemberAt(co,ip,el) - subLi0 = constitutive_mech_partionedLi0(ph)%data(1:3,1:3,me) - subLp0 = crystallite_partitionedLp0(1:3,1:3,co,ip,el) - plasticState (material_phaseAt(co,el))%subState0( :,material_phaseMemberAt(co,ip,el)) = & - plasticState (material_phaseAt(co,el))%partitionedState0(:,material_phaseMemberAt(co,ip,el)) - - do s = 1, phase_Nsources(material_phaseAt(co,el)) - sourceState(material_phaseAt(co,el))%p(s)%subState0( :,material_phaseMemberAt(co,ip,el)) = & - sourceState(material_phaseAt(co,el))%p(s)%partitionedState0(:,material_phaseMemberAt(co,ip,el)) - enddo - crystallite_subFp0(1:3,1:3,co,ip,el) = constitutive_mech_partionedFp0(ph)%data(1:3,1:3,me) - crystallite_subFi0(1:3,1:3,co,ip,el) = constitutive_mech_partionedFi0(ph)%data(1:3,1:3,me) - crystallite_subF0(1:3,1:3,co,ip,el) = crystallite_partitionedF0(1:3,1:3,co,ip,el) - subFrac = 0.0_pReal - subStep = 1.0_pReal/num%subStepSizeCryst - todo = .true. - crystallite_converged(co,ip,el) = .false. ! pretend failed step of 1/subStepSizeCryst - - todo = .true. - NiterationCrystallite = 0 - cutbackLooping: do while (todo) - NiterationCrystallite = NiterationCrystallite + 1 - -!-------------------------------------------------------------------------------------------------- -! wind forward - if (crystallite_converged(co,ip,el)) then - formerSubStep = subStep - subFrac = subFrac + subStep - subStep = min(1.0_pReal - subFrac, num%stepIncreaseCryst * subStep) - - todo = subStep > 0.0_pReal ! still time left to integrate on? - - if (todo) then - crystallite_subF0 (1:3,1:3,co,ip,el) = crystallite_subF(1:3,1:3,co,ip,el) - subLp0 = crystallite_Lp (1:3,1:3,co,ip,el) - subLi0 = constitutive_mech_Li(ph)%data(1:3,1:3,me) - crystallite_subFp0(1:3,1:3,co,ip,el) = constitutive_mech_Fp(ph)%data(1:3,1:3,me) - crystallite_subFi0(1:3,1:3,co,ip,el) = constitutive_mech_Fi(ph)%data(1:3,1:3,me) - plasticState( material_phaseAt(co,el))%subState0(:,material_phaseMemberAt(co,ip,el)) & - = plasticState(material_phaseAt(co,el))%state( :,material_phaseMemberAt(co,ip,el)) - do s = 1, phase_Nsources(material_phaseAt(co,el)) - sourceState( material_phaseAt(co,el))%p(s)%subState0(:,material_phaseMemberAt(co,ip,el)) & - = sourceState(material_phaseAt(co,el))%p(s)%state( :,material_phaseMemberAt(co,ip,el)) - enddo - endif -!-------------------------------------------------------------------------------------------------- -! cut back (reduced time and restore) - else - subStep = num%subStepSizeCryst * subStep - constitutive_mech_Fp(ph)%data(1:3,1:3,me) = crystallite_subFp0(1:3,1:3,co,ip,el) - constitutive_mech_Fi(ph)%data(1:3,1:3,me) = crystallite_subFi0(1:3,1:3,co,ip,el) - crystallite_S (1:3,1:3,co,ip,el) = crystallite_S0 (1:3,1:3,co,ip,el) - if (subStep < 1.0_pReal) then ! actual (not initial) cutback - crystallite_Lp (1:3,1:3,co,ip,el) = subLp0 - constitutive_mech_Li(ph)%data(1:3,1:3,me) = subLi0 - endif - plasticState (material_phaseAt(co,el))%state( :,material_phaseMemberAt(co,ip,el)) & - = plasticState(material_phaseAt(co,el))%subState0(:,material_phaseMemberAt(co,ip,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)%subState0(:,material_phaseMemberAt(co,ip,el)) - enddo - - ! cant restore dotState here, since not yet calculated in first cutback after initialization - todo = subStep > num%subStepMinCryst ! still on track or already done (beyond repair) - endif - -!-------------------------------------------------------------------------------------------------- -! prepare for integration - if (todo) then - crystallite_subF(1:3,1:3,co,ip,el) = crystallite_subF0(1:3,1:3,co,ip,el) & - + subStep *( crystallite_partitionedF (1:3,1:3,co,ip,el) & - -crystallite_partitionedF0(1:3,1:3,co,ip,el)) - crystallite_Fe(1:3,1:3,co,ip,el) = matmul(crystallite_subF(1:3,1:3,co,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(co,ip,el) = subStep * dt - crystallite_converged(co,ip,el) = .false. - call integrateState(co,ip,el) - call integrateSourceState(co,ip,el) - endif - - enddo cutbackLooping - -! return whether converged or not - crystallite_stress = crystallite_converged(co,ip,el) - -end function crystallite_stress !-------------------------------------------------------------------------------------------------- diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90 index 6a60b4ade..d448b0505 100644 --- a/src/constitutive_mech.f90 +++ b/src/constitutive_mech.f90 @@ -290,6 +290,8 @@ submodule(constitutive) constitutive_mech end type tOutput type(tOutput), allocatable, dimension(:) :: output_constituent + procedure(integrateStateFPI), pointer :: integrateState + contains @@ -959,7 +961,7 @@ 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(co,ip,el) +subroutine integrateStateFPI(co,ip,el) integer, intent(in) :: & el, & !< element index in element loop @@ -1475,5 +1477,121 @@ module function constitutive_homogenizedC(co,ip,el) result(C) end function constitutive_homogenizedC + +!-------------------------------------------------------------------------------------------------- +!> @brief calculate stress (P) +!-------------------------------------------------------------------------------------------------- +module function crystallite_stress(dt,co,ip,el) + + real(pReal), intent(in) :: dt + integer, intent(in) :: & + co, & + ip, & + el + + logical :: crystallite_stress + + real(pReal) :: & + formerSubStep + integer :: & + NiterationCrystallite, & ! number of iterations in crystallite loop + s, ph, me + logical :: todo + real(pReal) :: subFrac,subStep + real(pReal), dimension(3,3) :: & + subLp0, & !< plastic velocity grad at start of crystallite inc + subLi0 !< intermediate velocity grad at start of crystallite inc + + + ph = material_phaseAt(co,el) + me = material_phaseMemberAt(co,ip,el) + subLi0 = constitutive_mech_partionedLi0(ph)%data(1:3,1:3,me) + subLp0 = crystallite_partitionedLp0(1:3,1:3,co,ip,el) + plasticState (material_phaseAt(co,el))%subState0( :,material_phaseMemberAt(co,ip,el)) = & + plasticState (material_phaseAt(co,el))%partitionedState0(:,material_phaseMemberAt(co,ip,el)) + + do s = 1, phase_Nsources(material_phaseAt(co,el)) + sourceState(material_phaseAt(co,el))%p(s)%subState0( :,material_phaseMemberAt(co,ip,el)) = & + sourceState(material_phaseAt(co,el))%p(s)%partitionedState0(:,material_phaseMemberAt(co,ip,el)) + enddo + crystallite_subFp0(1:3,1:3,co,ip,el) = constitutive_mech_partionedFp0(ph)%data(1:3,1:3,me) + crystallite_subFi0(1:3,1:3,co,ip,el) = constitutive_mech_partionedFi0(ph)%data(1:3,1:3,me) + crystallite_subF0(1:3,1:3,co,ip,el) = crystallite_partitionedF0(1:3,1:3,co,ip,el) + subFrac = 0.0_pReal + subStep = 1.0_pReal/num%subStepSizeCryst + todo = .true. + crystallite_converged(co,ip,el) = .false. ! pretend failed step of 1/subStepSizeCryst + + todo = .true. + NiterationCrystallite = 0 + cutbackLooping: do while (todo) + NiterationCrystallite = NiterationCrystallite + 1 + +!-------------------------------------------------------------------------------------------------- +! wind forward + if (crystallite_converged(co,ip,el)) then + formerSubStep = subStep + subFrac = subFrac + subStep + subStep = min(1.0_pReal - subFrac, num%stepIncreaseCryst * subStep) + + todo = subStep > 0.0_pReal ! still time left to integrate on? + + if (todo) then + crystallite_subF0 (1:3,1:3,co,ip,el) = crystallite_subF(1:3,1:3,co,ip,el) + subLp0 = crystallite_Lp (1:3,1:3,co,ip,el) + subLi0 = constitutive_mech_Li(ph)%data(1:3,1:3,me) + crystallite_subFp0(1:3,1:3,co,ip,el) = constitutive_mech_Fp(ph)%data(1:3,1:3,me) + crystallite_subFi0(1:3,1:3,co,ip,el) = constitutive_mech_Fi(ph)%data(1:3,1:3,me) + plasticState( material_phaseAt(co,el))%subState0(:,material_phaseMemberAt(co,ip,el)) & + = plasticState(material_phaseAt(co,el))%state( :,material_phaseMemberAt(co,ip,el)) + do s = 1, phase_Nsources(material_phaseAt(co,el)) + sourceState( material_phaseAt(co,el))%p(s)%subState0(:,material_phaseMemberAt(co,ip,el)) & + = sourceState(material_phaseAt(co,el))%p(s)%state( :,material_phaseMemberAt(co,ip,el)) + enddo + endif +!-------------------------------------------------------------------------------------------------- +! cut back (reduced time and restore) + else + subStep = num%subStepSizeCryst * subStep + constitutive_mech_Fp(ph)%data(1:3,1:3,me) = crystallite_subFp0(1:3,1:3,co,ip,el) + constitutive_mech_Fi(ph)%data(1:3,1:3,me) = crystallite_subFi0(1:3,1:3,co,ip,el) + crystallite_S (1:3,1:3,co,ip,el) = crystallite_S0 (1:3,1:3,co,ip,el) + if (subStep < 1.0_pReal) then ! actual (not initial) cutback + crystallite_Lp (1:3,1:3,co,ip,el) = subLp0 + constitutive_mech_Li(ph)%data(1:3,1:3,me) = subLi0 + endif + plasticState (material_phaseAt(co,el))%state( :,material_phaseMemberAt(co,ip,el)) & + = plasticState(material_phaseAt(co,el))%subState0(:,material_phaseMemberAt(co,ip,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)%subState0(:,material_phaseMemberAt(co,ip,el)) + enddo + + ! cant restore dotState here, since not yet calculated in first cutback after initialization + todo = subStep > num%subStepMinCryst ! still on track or already done (beyond repair) + endif + +!-------------------------------------------------------------------------------------------------- +! prepare for integration + if (todo) then + crystallite_subF(1:3,1:3,co,ip,el) = crystallite_subF0(1:3,1:3,co,ip,el) & + + subStep *( crystallite_partitionedF (1:3,1:3,co,ip,el) & + -crystallite_partitionedF0(1:3,1:3,co,ip,el)) + crystallite_Fe(1:3,1:3,co,ip,el) = matmul(crystallite_subF(1:3,1:3,co,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(co,ip,el) = subStep * dt + crystallite_converged(co,ip,el) = .false. + call integrateState(co,ip,el) + call integrateSourceState(co,ip,el) + endif + + enddo cutbackLooping + +! return whether converged or not + crystallite_stress = crystallite_converged(co,ip,el) + +end function crystallite_stress + end submodule constitutive_mech