From 5fce37fb3e4effdc97cdac271add21896a912550 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 21 Dec 2020 10:57:18 +0100 Subject: [PATCH] only relevant for mechanics --- src/constitutive.f90 | 600 ++------------------------------------ src/constitutive_mech.f90 | 576 ++++++++++++++++++++++++++++++++++++ 2 files changed, 599 insertions(+), 577 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 75d5e098f..58442283a 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -392,6 +392,26 @@ end function constitutive_deltaState 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 + end subroutine integrateStateFPI + + module subroutine integrateStateEuler(g,i,e) + integer, intent(in) :: e, i, g + end subroutine integrateStateEuler + + module subroutine integrateStateAdaptiveEuler(g,i,e) + integer, intent(in) :: e, i, g + end subroutine integrateStateAdaptiveEuler + + module subroutine integrateStateRK4(g,i,e) + integer, intent(in) :: e, i, g + end subroutine integrateStateRK4 + + module subroutine integrateStateRKCK45(g,i,e) + integer, intent(in) :: e, i, g + end subroutine integrateStateRKCK45 + end interface @@ -414,9 +434,8 @@ end function constitutive_deltaState plastic_nonlocal_updateCompatibility, & plastic_active, & source_active, & - kinematics_active - - public :: & + kinematics_active, & + converged, & crystallite_init, & crystallite_stress, & crystallite_stressTangent, & @@ -429,6 +448,7 @@ end function constitutive_deltaState crystallite_initializeRestorationPoints, & crystallite_windForward, & crystallite_restore + contains @@ -1562,338 +1582,6 @@ subroutine crystallite_results end subroutine crystallite_results -!-------------------------------------------------------------------------------------------------- -!> @brief calculation of stress (P) with time integration based on a residuum in Lp and -!> intermediate acceleration of the Newton-Raphson correction -!-------------------------------------------------------------------------------------------------- -function integrateStress(ipc,ip,el,timeFraction) result(broken) - - integer, intent(in):: el, & ! element index - ip, & ! integration point index - ipc ! grain index - real(pReal), optional, intent(in) :: timeFraction ! fraction of timestep - - real(pReal), dimension(3,3):: F, & ! deformation gradient at end of timestep - Fp_new, & ! plastic deformation gradient at end of timestep - invFp_new, & ! inverse of Fp_new - invFp_current, & ! inverse of Fp_current - Lpguess, & ! current guess for plastic velocity gradient - Lpguess_old, & ! known last good guess for plastic velocity gradient - Lp_constitutive, & ! plastic velocity gradient resulting from constitutive law - residuumLp, & ! current residuum of plastic velocity gradient - residuumLp_old, & ! last residuum of plastic velocity gradient - deltaLp, & ! direction of next guess - Fi_new, & ! gradient of intermediate deformation stages - invFi_new, & - invFi_current, & ! inverse of Fi_current - Liguess, & ! current guess for intermediate velocity gradient - Liguess_old, & ! known last good guess for intermediate velocity gradient - Li_constitutive, & ! intermediate velocity gradient resulting from constitutive law - residuumLi, & ! current residuum of intermediate velocity gradient - residuumLi_old, & ! last residuum of intermediate velocity gradient - deltaLi, & ! direction of next guess - Fe, & ! elastic deformation gradient - S, & ! 2nd Piola-Kirchhoff Stress in plastic (lattice) configuration - A, & - B, & - temp_33 - real(pReal), dimension(9) :: temp_9 ! needed for matrix inversion by LAPACK - integer, dimension(9) :: devNull_9 ! needed for matrix inversion by LAPACK - real(pReal), dimension(9,9) :: dRLp_dLp, & ! partial derivative of residuum (Jacobian for Newton-Raphson scheme) - dRLi_dLi ! partial derivative of residuumI (Jacobian for Newton-Raphson scheme) - real(pReal), dimension(3,3,3,3):: dS_dFe, & ! partial derivative of 2nd Piola-Kirchhoff stress - dS_dFi, & - dFe_dLp, & ! partial derivative of elastic deformation gradient - dFe_dLi, & - dFi_dLi, & - dLp_dFi, & - dLi_dFi, & - dLp_dS, & - dLi_dS - real(pReal) steplengthLp, & - steplengthLi, & - dt, & ! time increment - atol_Lp, & - atol_Li, & - devNull - integer NiterationStressLp, & ! number of stress integrations - NiterationStressLi, & ! number of inner stress integrations - ierr, & ! error indicator for LAPACK - o, & - p, & - m, & - jacoCounterLp, & - jacoCounterLi ! counters to check for Jacobian update - logical :: error,broken - - broken = .true. - - if (present(timeFraction)) then - dt = crystallite_subdt(ipc,ip,el) * timeFraction - F = crystallite_subF0(1:3,1:3,ipc,ip,el) & - + (crystallite_subF(1:3,1:3,ipc,ip,el) - crystallite_subF0(1:3,1:3,ipc,ip,el)) * timeFraction - else - dt = crystallite_subdt(ipc,ip,el) - F = crystallite_subF(1:3,1:3,ipc,ip,el) - endif - - call constitutive_plastic_dependentState(crystallite_partitionedF(1:3,1:3,ipc,ip,el), & - crystallite_Fp(1:3,1:3,ipc,ip,el),ipc,ip,el) - - p = material_phaseAt(ipc,el) - m = material_phaseMemberAt(ipc,ip,el) - - Lpguess = crystallite_Lp(1:3,1:3,ipc,ip,el) ! take as first guess - Liguess = constitutive_mech_Li(p)%data(1:3,1:3,m) ! take as first guess - - call math_invert33(invFp_current,devNull,error,crystallite_subFp0(1:3,1:3,ipc,ip,el)) - if (error) return ! error - call math_invert33(invFi_current,devNull,error,crystallite_subFi0(1:3,1:3,ipc,ip,el)) - if (error) return ! error - - A = matmul(F,invFp_current) ! intermediate tensor needed later to calculate dFe_dLp - - jacoCounterLi = 0 - steplengthLi = 1.0_pReal - residuumLi_old = 0.0_pReal - Liguess_old = Liguess - - NiterationStressLi = 0 - LiLoop: do - NiterationStressLi = NiterationStressLi + 1 - if (NiterationStressLi>num%nStress) return ! error - - invFi_new = matmul(invFi_current,math_I3 - dt*Liguess) - Fi_new = math_inv33(invFi_new) - - jacoCounterLp = 0 - steplengthLp = 1.0_pReal - residuumLp_old = 0.0_pReal - Lpguess_old = Lpguess - - NiterationStressLp = 0 - LpLoop: do - NiterationStressLp = NiterationStressLp + 1 - if (NiterationStressLp>num%nStress) return ! error - - B = math_I3 - dt*Lpguess - Fe = matmul(matmul(A,B), invFi_new) - call constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & - Fe, Fi_new, ipc, ip, el) - - call constitutive_plastic_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, & - S, Fi_new, ipc, ip, el) - - !* update current residuum and check for convergence of loop - atol_Lp = max(num%rtol_crystalliteStress * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error - num%atol_crystalliteStress) ! minimum lower cutoff - residuumLp = Lpguess - Lp_constitutive - - if (any(IEEE_is_NaN(residuumLp))) then - return ! error - elseif (norm2(residuumLp) < atol_Lp) then ! converged if below absolute tolerance - exit LpLoop - elseif (NiterationStressLp == 1 .or. norm2(residuumLp) < norm2(residuumLp_old)) then ! not converged, but improved norm of residuum (always proceed in first iteration)... - residuumLp_old = residuumLp ! ...remember old values and... - Lpguess_old = Lpguess - steplengthLp = 1.0_pReal ! ...proceed with normal step length (calculate new search direction) - else ! not converged and residuum not improved... - steplengthLp = num%subStepSizeLp * steplengthLp ! ...try with smaller step length in same direction - Lpguess = Lpguess_old & - + deltaLp * stepLengthLp - cycle LpLoop - endif - - calculateJacobiLi: if (mod(jacoCounterLp, num%iJacoLpresiduum) == 0) then - jacoCounterLp = jacoCounterLp + 1 - - do o=1,3; do p=1,3 - dFe_dLp(o,1:3,p,1:3) = - dt * A(o,p)*transpose(invFi_new) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j) - enddo; enddo - dRLp_dLp = math_eye(9) & - - math_3333to99(math_mul3333xx3333(math_mul3333xx3333(dLp_dS,dS_dFe),dFe_dLp)) - temp_9 = math_33to9(residuumLp) - call dgesv(9,1,dRLp_dLp,9,devNull_9,temp_9,9,ierr) ! solve dRLp/dLp * delta Lp = -res for delta Lp - if (ierr /= 0) return ! error - deltaLp = - math_9to33(temp_9) - endif calculateJacobiLi - - Lpguess = Lpguess & - + deltaLp * steplengthLp - enddo LpLoop - - call constitutive_LiAndItsTangents(Li_constitutive, dLi_dS, dLi_dFi, & - S, Fi_new, ipc, ip, el) - - !* update current residuum and check for convergence of loop - atol_Li = max(num%rtol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error - num%atol_crystalliteStress) ! minimum lower cutoff - residuumLi = Liguess - Li_constitutive - if (any(IEEE_is_NaN(residuumLi))) then - return ! error - elseif (norm2(residuumLi) < atol_Li) then ! converged if below absolute tolerance - exit LiLoop - elseif (NiterationStressLi == 1 .or. norm2(residuumLi) < norm2(residuumLi_old)) then ! not converged, but improved norm of residuum (always proceed in first iteration)... - residuumLi_old = residuumLi ! ...remember old values and... - Liguess_old = Liguess - steplengthLi = 1.0_pReal ! ...proceed with normal step length (calculate new search direction) - else ! not converged and residuum not improved... - steplengthLi = num%subStepSizeLi * steplengthLi ! ...try with smaller step length in same direction - Liguess = Liguess_old & - + deltaLi * steplengthLi - cycle LiLoop - endif - - calculateJacobiLp: if (mod(jacoCounterLi, num%iJacoLpresiduum) == 0) then - jacoCounterLi = jacoCounterLi + 1 - - temp_33 = matmul(matmul(A,B),invFi_current) - do o=1,3; do p=1,3 - dFe_dLi(1:3,o,1:3,p) = -dt*math_I3(o,p)*temp_33 ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j) - dFi_dLi(1:3,o,1:3,p) = -dt*math_I3(o,p)*invFi_current - enddo; enddo - do o=1,3; do p=1,3 - dFi_dLi(1:3,1:3,o,p) = matmul(matmul(Fi_new,dFi_dLi(1:3,1:3,o,p)),Fi_new) - enddo; enddo - dRLi_dLi = math_eye(9) & - - math_3333to99(math_mul3333xx3333(dLi_dS, math_mul3333xx3333(dS_dFe, dFe_dLi) & - + math_mul3333xx3333(dS_dFi, dFi_dLi))) & - - math_3333to99(math_mul3333xx3333(dLi_dFi, dFi_dLi)) - temp_9 = math_33to9(residuumLi) - call dgesv(9,1,dRLi_dLi,9,devNull_9,temp_9,9,ierr) ! solve dRLi/dLp * delta Li = -res for delta Li - if (ierr /= 0) return ! error - deltaLi = - math_9to33(temp_9) - endif calculateJacobiLp - - Liguess = Liguess & - + deltaLi * steplengthLi - enddo LiLoop - - invFp_new = matmul(invFp_current,B) - call math_invert33(Fp_new,devNull,error,invFp_new) - if (error) return ! error - - p = material_phaseAt(ipc,el) - m = material_phaseMemberAt(ipc,ip,el) - - crystallite_P (1:3,1:3,ipc,ip,el) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new))) - crystallite_S (1:3,1:3,ipc,ip,el) = S - crystallite_Lp (1:3,1:3,ipc,ip,el) = Lpguess - constitutive_mech_Li(p)%data(1:3,1:3,m) = Liguess - crystallite_Fp (1:3,1:3,ipc,ip,el) = Fp_new / math_det33(Fp_new)**(1.0_pReal/3.0_pReal) ! regularize - constitutive_mech_Fi(p)%data(1:3,1:3,m) = Fi_new - crystallite_Fe (1:3,1:3,ipc,ip,el) = matmul(matmul(F,invFp_new),invFi_new) - broken = .false. - -end function integrateStress - - -!-------------------------------------------------------------------------------------------------- -!> @brief integrate stress, state with adaptive 1st order explicit Euler method -!> using Fixed Point Iteration to adapt the stepsize -!-------------------------------------------------------------------------------------------------- -subroutine integrateStateFPI(g,i,e) - - integer, intent(in) :: & - e, & !< element index in element loop - i, & !< integration point index in ip loop - g !< grain index in grain loop - integer :: & - NiterationState, & !< number of iterations in state loop - p, & - c, & - s, & - size_pl - integer, dimension(maxval(phase_Nsources)) :: & - size_so - real(pReal) :: & - zeta - real(pReal), dimension(max(constitutive_plasticity_maxSizeDotState,constitutive_source_maxSizeDotState)) :: & - r ! state residuum - real(pReal), dimension(constitutive_plasticity_maxSizeDotState,2) :: & - plastic_dotState - real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState - logical :: & - broken - - p = material_phaseAt(g,e) - c = material_phaseMemberAt(g,i,e) - - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partitionedF0, & - constitutive_mech_Fi(p)%data(1:3,1:3,c), & - crystallite_partitionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) - 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) - 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) - - broken = integrateStress(g,i,e) - if(broken) exit iteration - - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partitionedF0, & - constitutive_mech_Fi(p)%data(1:3,1:3,c), & - crystallite_partitionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken) exit iteration - - zeta = damper(plasticState(p)%dotState(:,c),plastic_dotState(1:size_pl,1),& - plastic_dotState(1:size_pl,2)) - plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * 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) - crystallite_converged(g,i,e) = converged(r(1:size_pl), & - plasticState(p)%state(1:size_pl,c), & - plasticState(p)%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) - exit iteration - endif - - enddo iteration - - - contains - - !-------------------------------------------------------------------------------------------------- - !> @brief calculate the damping for correction of state and dot state - !-------------------------------------------------------------------------------------------------- - real(pReal) pure function damper(current,previous,previous2) - - real(pReal), dimension(:), intent(in) ::& - current, previous, previous2 - - real(pReal) :: dot_prod12, dot_prod22 - - dot_prod12 = dot_product(current - previous, previous - previous2) - dot_prod22 = dot_product(previous - previous2, previous - previous2) - if ((dot_product(current,previous) < 0.0_pReal .or. dot_prod12 < 0.0_pReal) .and. dot_prod22 > 0.0_pReal) then - damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22) - else - damper = 1.0_pReal - endif - - end function damper - -end subroutine integrateStateFPI - - !-------------------------------------------------------------------------------------------------- !> @brief integrate stress, state with adaptive 1st order explicit Euler method !> using Fixed Point Iteration to adapt the stepsize @@ -1993,248 +1681,6 @@ subroutine integrateSourceState(g,i,e) end subroutine integrateSourceState -!-------------------------------------------------------------------------------------------------- -!> @brief integrate state with 1st order explicit Euler method -!-------------------------------------------------------------------------------------------------- -subroutine integrateStateEuler(g,i,e) - - integer, intent(in) :: & - e, & !< element index in element loop - i, & !< integration point index in ip loop - g !< grain index in grain loop - integer :: & - p, & - c, & - sizeDotState - logical :: & - broken - - p = material_phaseAt(g,e) - c = material_phaseMemberAt(g,i,e) - - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partitionedF0, & - constitutive_mech_Fi(p)%data(1:3,1:3,c), & - crystallite_partitionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) - 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) - - 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(broken) return - - broken = integrateStress(g,i,e) - crystallite_converged(g,i,e) = .not. broken - -end subroutine integrateStateEuler - - -!-------------------------------------------------------------------------------------------------- -!> @brief integrate stress, state with 1st order Euler method with adaptive step size -!-------------------------------------------------------------------------------------------------- -subroutine integrateStateAdaptiveEuler(g,i,e) - - integer, intent(in) :: & - e, & !< element index in element loop - i, & !< integration point index in ip loop - g !< grain index in grain loop - integer :: & - p, & - c, & - sizeDotState - logical :: & - broken - - real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: residuum_plastic - - - p = material_phaseAt(g,e) - c = material_phaseMemberAt(g,i,e) - - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partitionedF0, & - constitutive_mech_Fi(p)%data(1:3,1:3,c), & - crystallite_partitionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken) return - - sizeDotState = plasticState(p)%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) - - 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(broken) return - - broken = integrateStress(g,i,e) - if(broken) return - - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partitionedF0, & - constitutive_mech_Fi(p)%data(1:3,1:3,c), & - crystallite_partitionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) - 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)) - -end subroutine integrateStateAdaptiveEuler - - -!--------------------------------------------------------------------------------------------------- -!> @brief Integrate state (including stress integration) with the classic Runge Kutta method -!--------------------------------------------------------------------------------------------------- -subroutine integrateStateRK4(g,i,e) - - integer, intent(in) :: g,i,e - - real(pReal), dimension(3,3), parameter :: & - A = reshape([& - 0.5_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.5_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, 1.0_pReal],& - shape(A)) - real(pReal), dimension(3), parameter :: & - C = [0.5_pReal, 0.5_pReal, 1.0_pReal] - 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) - -end subroutine integrateStateRK4 - - -!--------------------------------------------------------------------------------------------------- -!> @brief Integrate state (including stress integration) with the Cash-Carp method -!--------------------------------------------------------------------------------------------------- -subroutine integrateStateRKCK45(g,i,e) - - integer, intent(in) :: g,i,e - - real(pReal), dimension(5,5), parameter :: & - A = reshape([& - 1._pReal/5._pReal, .0_pReal, .0_pReal, .0_pReal, .0_pReal, & - 3._pReal/40._pReal, 9._pReal/40._pReal, .0_pReal, .0_pReal, .0_pReal, & - 3_pReal/10._pReal, -9._pReal/10._pReal, 6._pReal/5._pReal, .0_pReal, .0_pReal, & - -11._pReal/54._pReal, 5._pReal/2._pReal, -70.0_pReal/27.0_pReal, 35.0_pReal/27.0_pReal, .0_pReal, & - 1631._pReal/55296._pReal,175._pReal/512._pReal,575._pReal/13824._pReal,44275._pReal/110592._pReal,253._pReal/4096._pReal],& - shape(A)) - real(pReal), dimension(5), parameter :: & - C = [0.2_pReal, 0.3_pReal, 0.6_pReal, 1.0_pReal, 0.875_pReal] - real(pReal), dimension(6), parameter :: & - B = & - [37.0_pReal/378.0_pReal, .0_pReal, 250.0_pReal/621.0_pReal, & - 125.0_pReal/594.0_pReal, .0_pReal, 512.0_pReal/1771.0_pReal], & - DB = B - & - [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) - -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) - - - real(pReal), dimension(:,:), intent(in) :: A - real(pReal), dimension(:), intent(in) :: B, CC - 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 - integer :: & - stage, & ! stage index in integration stage loop - n, & - p, & - c, & - sizeDotState - logical :: & - broken - real(pReal), dimension(constitutive_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState - - p = material_phaseAt(g,e) - c = material_phaseMemberAt(g,i,e) - - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partitionedF0, & - constitutive_mech_Fi(p)%data(1:3,1:3,c), & - crystallite_partitionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) - 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) - - do n = 2, stage - sizeDotState = plasticState(p)%sizeDotState - plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) & - + 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) - - broken = integrateStress(g,i,e,CC(stage)) - if(broken) exit - - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partitionedF0, & - constitutive_mech_Fi(p)%data(1:3,1:3,c), & - crystallite_partitionedFp0, & - crystallite_subdt(g,i,e)*CC(stage), g,i,e,p,c) - if(broken) exit - - enddo - if(broken) return - - sizeDotState = plasticState(p)%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) - 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)) - - 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) - if(broken) return - - broken = integrateStress(g,i,e) - crystallite_converged(g,i,e) = .not. broken - - -end subroutine integrateStateRK - !-------------------------------------------------------------------------------------------------- !> @brief determines whether a point is converged diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90 index 8f08aa08e..dea5ed647 100644 --- a/src/constitutive_mech.f90 +++ b/src/constitutive_mech.f90 @@ -691,5 +691,581 @@ module subroutine plastic_results end subroutine plastic_results + +!-------------------------------------------------------------------------------------------------- +!> @brief calculation of stress (P) with time integration based on a residuum in Lp and +!> intermediate acceleration of the Newton-Raphson correction +!-------------------------------------------------------------------------------------------------- +function integrateStress(ipc,ip,el,timeFraction) result(broken) + + integer, intent(in):: el, & ! element index + ip, & ! integration point index + ipc ! grain index + real(pReal), optional, intent(in) :: timeFraction ! fraction of timestep + + real(pReal), dimension(3,3):: F, & ! deformation gradient at end of timestep + Fp_new, & ! plastic deformation gradient at end of timestep + invFp_new, & ! inverse of Fp_new + invFp_current, & ! inverse of Fp_current + Lpguess, & ! current guess for plastic velocity gradient + Lpguess_old, & ! known last good guess for plastic velocity gradient + Lp_constitutive, & ! plastic velocity gradient resulting from constitutive law + residuumLp, & ! current residuum of plastic velocity gradient + residuumLp_old, & ! last residuum of plastic velocity gradient + deltaLp, & ! direction of next guess + Fi_new, & ! gradient of intermediate deformation stages + invFi_new, & + invFi_current, & ! inverse of Fi_current + Liguess, & ! current guess for intermediate velocity gradient + Liguess_old, & ! known last good guess for intermediate velocity gradient + Li_constitutive, & ! intermediate velocity gradient resulting from constitutive law + residuumLi, & ! current residuum of intermediate velocity gradient + residuumLi_old, & ! last residuum of intermediate velocity gradient + deltaLi, & ! direction of next guess + Fe, & ! elastic deformation gradient + S, & ! 2nd Piola-Kirchhoff Stress in plastic (lattice) configuration + A, & + B, & + temp_33 + real(pReal), dimension(9) :: temp_9 ! needed for matrix inversion by LAPACK + integer, dimension(9) :: devNull_9 ! needed for matrix inversion by LAPACK + real(pReal), dimension(9,9) :: dRLp_dLp, & ! partial derivative of residuum (Jacobian for Newton-Raphson scheme) + dRLi_dLi ! partial derivative of residuumI (Jacobian for Newton-Raphson scheme) + real(pReal), dimension(3,3,3,3):: dS_dFe, & ! partial derivative of 2nd Piola-Kirchhoff stress + dS_dFi, & + dFe_dLp, & ! partial derivative of elastic deformation gradient + dFe_dLi, & + dFi_dLi, & + dLp_dFi, & + dLi_dFi, & + dLp_dS, & + dLi_dS + real(pReal) steplengthLp, & + steplengthLi, & + dt, & ! time increment + atol_Lp, & + atol_Li, & + devNull + integer NiterationStressLp, & ! number of stress integrations + NiterationStressLi, & ! number of inner stress integrations + ierr, & ! error indicator for LAPACK + o, & + p, & + m, & + jacoCounterLp, & + jacoCounterLi ! counters to check for Jacobian update + logical :: error,broken + + broken = .true. + + if (present(timeFraction)) then + dt = crystallite_subdt(ipc,ip,el) * timeFraction + F = crystallite_subF0(1:3,1:3,ipc,ip,el) & + + (crystallite_subF(1:3,1:3,ipc,ip,el) - crystallite_subF0(1:3,1:3,ipc,ip,el)) * timeFraction + else + dt = crystallite_subdt(ipc,ip,el) + F = crystallite_subF(1:3,1:3,ipc,ip,el) + endif + + call constitutive_plastic_dependentState(crystallite_partitionedF(1:3,1:3,ipc,ip,el), & + crystallite_Fp(1:3,1:3,ipc,ip,el),ipc,ip,el) + + p = material_phaseAt(ipc,el) + m = material_phaseMemberAt(ipc,ip,el) + + Lpguess = crystallite_Lp(1:3,1:3,ipc,ip,el) ! take as first guess + Liguess = constitutive_mech_Li(p)%data(1:3,1:3,m) ! take as first guess + + call math_invert33(invFp_current,devNull,error,crystallite_subFp0(1:3,1:3,ipc,ip,el)) + if (error) return ! error + call math_invert33(invFi_current,devNull,error,crystallite_subFi0(1:3,1:3,ipc,ip,el)) + if (error) return ! error + + A = matmul(F,invFp_current) ! intermediate tensor needed later to calculate dFe_dLp + + jacoCounterLi = 0 + steplengthLi = 1.0_pReal + residuumLi_old = 0.0_pReal + Liguess_old = Liguess + + NiterationStressLi = 0 + LiLoop: do + NiterationStressLi = NiterationStressLi + 1 + if (NiterationStressLi>num%nStress) return ! error + + invFi_new = matmul(invFi_current,math_I3 - dt*Liguess) + Fi_new = math_inv33(invFi_new) + + jacoCounterLp = 0 + steplengthLp = 1.0_pReal + residuumLp_old = 0.0_pReal + Lpguess_old = Lpguess + + NiterationStressLp = 0 + LpLoop: do + NiterationStressLp = NiterationStressLp + 1 + if (NiterationStressLp>num%nStress) return ! error + + B = math_I3 - dt*Lpguess + Fe = matmul(matmul(A,B), invFi_new) + call constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & + Fe, Fi_new, ipc, ip, el) + + call constitutive_plastic_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, & + S, Fi_new, ipc, ip, el) + + !* update current residuum and check for convergence of loop + atol_Lp = max(num%rtol_crystalliteStress * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error + num%atol_crystalliteStress) ! minimum lower cutoff + residuumLp = Lpguess - Lp_constitutive + + if (any(IEEE_is_NaN(residuumLp))) then + return ! error + elseif (norm2(residuumLp) < atol_Lp) then ! converged if below absolute tolerance + exit LpLoop + elseif (NiterationStressLp == 1 .or. norm2(residuumLp) < norm2(residuumLp_old)) then ! not converged, but improved norm of residuum (always proceed in first iteration)... + residuumLp_old = residuumLp ! ...remember old values and... + Lpguess_old = Lpguess + steplengthLp = 1.0_pReal ! ...proceed with normal step length (calculate new search direction) + else ! not converged and residuum not improved... + steplengthLp = num%subStepSizeLp * steplengthLp ! ...try with smaller step length in same direction + Lpguess = Lpguess_old & + + deltaLp * stepLengthLp + cycle LpLoop + endif + + calculateJacobiLi: if (mod(jacoCounterLp, num%iJacoLpresiduum) == 0) then + jacoCounterLp = jacoCounterLp + 1 + + do o=1,3; do p=1,3 + dFe_dLp(o,1:3,p,1:3) = - dt * A(o,p)*transpose(invFi_new) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j) + enddo; enddo + dRLp_dLp = math_eye(9) & + - math_3333to99(math_mul3333xx3333(math_mul3333xx3333(dLp_dS,dS_dFe),dFe_dLp)) + temp_9 = math_33to9(residuumLp) + call dgesv(9,1,dRLp_dLp,9,devNull_9,temp_9,9,ierr) ! solve dRLp/dLp * delta Lp = -res for delta Lp + if (ierr /= 0) return ! error + deltaLp = - math_9to33(temp_9) + endif calculateJacobiLi + + Lpguess = Lpguess & + + deltaLp * steplengthLp + enddo LpLoop + + call constitutive_LiAndItsTangents(Li_constitutive, dLi_dS, dLi_dFi, & + S, Fi_new, ipc, ip, el) + + !* update current residuum and check for convergence of loop + atol_Li = max(num%rtol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error + num%atol_crystalliteStress) ! minimum lower cutoff + residuumLi = Liguess - Li_constitutive + if (any(IEEE_is_NaN(residuumLi))) then + return ! error + elseif (norm2(residuumLi) < atol_Li) then ! converged if below absolute tolerance + exit LiLoop + elseif (NiterationStressLi == 1 .or. norm2(residuumLi) < norm2(residuumLi_old)) then ! not converged, but improved norm of residuum (always proceed in first iteration)... + residuumLi_old = residuumLi ! ...remember old values and... + Liguess_old = Liguess + steplengthLi = 1.0_pReal ! ...proceed with normal step length (calculate new search direction) + else ! not converged and residuum not improved... + steplengthLi = num%subStepSizeLi * steplengthLi ! ...try with smaller step length in same direction + Liguess = Liguess_old & + + deltaLi * steplengthLi + cycle LiLoop + endif + + calculateJacobiLp: if (mod(jacoCounterLi, num%iJacoLpresiduum) == 0) then + jacoCounterLi = jacoCounterLi + 1 + + temp_33 = matmul(matmul(A,B),invFi_current) + do o=1,3; do p=1,3 + dFe_dLi(1:3,o,1:3,p) = -dt*math_I3(o,p)*temp_33 ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j) + dFi_dLi(1:3,o,1:3,p) = -dt*math_I3(o,p)*invFi_current + enddo; enddo + do o=1,3; do p=1,3 + dFi_dLi(1:3,1:3,o,p) = matmul(matmul(Fi_new,dFi_dLi(1:3,1:3,o,p)),Fi_new) + enddo; enddo + dRLi_dLi = math_eye(9) & + - math_3333to99(math_mul3333xx3333(dLi_dS, math_mul3333xx3333(dS_dFe, dFe_dLi) & + + math_mul3333xx3333(dS_dFi, dFi_dLi))) & + - math_3333to99(math_mul3333xx3333(dLi_dFi, dFi_dLi)) + temp_9 = math_33to9(residuumLi) + call dgesv(9,1,dRLi_dLi,9,devNull_9,temp_9,9,ierr) ! solve dRLi/dLp * delta Li = -res for delta Li + if (ierr /= 0) return ! error + deltaLi = - math_9to33(temp_9) + endif calculateJacobiLp + + Liguess = Liguess & + + deltaLi * steplengthLi + enddo LiLoop + + invFp_new = matmul(invFp_current,B) + call math_invert33(Fp_new,devNull,error,invFp_new) + if (error) return ! error + + p = material_phaseAt(ipc,el) + m = material_phaseMemberAt(ipc,ip,el) + + crystallite_P (1:3,1:3,ipc,ip,el) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new))) + crystallite_S (1:3,1:3,ipc,ip,el) = S + crystallite_Lp (1:3,1:3,ipc,ip,el) = Lpguess + constitutive_mech_Li(p)%data(1:3,1:3,m) = Liguess + crystallite_Fp (1:3,1:3,ipc,ip,el) = Fp_new / math_det33(Fp_new)**(1.0_pReal/3.0_pReal) ! regularize + constitutive_mech_Fi(p)%data(1:3,1:3,m) = Fi_new + crystallite_Fe (1:3,1:3,ipc,ip,el) = matmul(matmul(F,invFp_new),invFi_new) + broken = .false. + +end function integrateStress + + +!-------------------------------------------------------------------------------------------------- +!> @brief integrate stress, state with adaptive 1st order explicit Euler method +!> using Fixed Point Iteration to adapt the stepsize +!-------------------------------------------------------------------------------------------------- +subroutine integrateStateFPI(g,i,e) + + integer, intent(in) :: & + e, & !< element index in element loop + i, & !< integration point index in ip loop + g !< grain index in grain loop + integer :: & + NiterationState, & !< number of iterations in state loop + p, & + c, & + s, & + size_pl + integer, dimension(maxval(phase_Nsources)) :: & + size_so + real(pReal) :: & + zeta + real(pReal), dimension(max(constitutive_plasticity_maxSizeDotState,constitutive_source_maxSizeDotState)) :: & + r ! state residuum + real(pReal), dimension(constitutive_plasticity_maxSizeDotState,2) :: & + plastic_dotState + real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState + logical :: & + broken + + p = material_phaseAt(g,e) + c = material_phaseMemberAt(g,i,e) + + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partitionedF0, & + constitutive_mech_Fi(p)%data(1:3,1:3,c), & + crystallite_partitionedFp0, & + crystallite_subdt(g,i,e), g,i,e,p,c) + 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) + 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) + + broken = integrateStress(g,i,e) + if(broken) exit iteration + + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partitionedF0, & + constitutive_mech_Fi(p)%data(1:3,1:3,c), & + crystallite_partitionedFp0, & + crystallite_subdt(g,i,e), g,i,e,p,c) + if(broken) exit iteration + + zeta = damper(plasticState(p)%dotState(:,c),plastic_dotState(1:size_pl,1),& + plastic_dotState(1:size_pl,2)) + plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * 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) + crystallite_converged(g,i,e) = converged(r(1:size_pl), & + plasticState(p)%state(1:size_pl,c), & + plasticState(p)%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) + exit iteration + endif + + enddo iteration + + + contains + + !-------------------------------------------------------------------------------------------------- + !> @brief calculate the damping for correction of state and dot state + !-------------------------------------------------------------------------------------------------- + real(pReal) pure function damper(current,previous,previous2) + + real(pReal), dimension(:), intent(in) ::& + current, previous, previous2 + + real(pReal) :: dot_prod12, dot_prod22 + + dot_prod12 = dot_product(current - previous, previous - previous2) + dot_prod22 = dot_product(previous - previous2, previous - previous2) + if ((dot_product(current,previous) < 0.0_pReal .or. dot_prod12 < 0.0_pReal) .and. dot_prod22 > 0.0_pReal) then + damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22) + else + damper = 1.0_pReal + endif + + end function damper + +end subroutine integrateStateFPI + + +!-------------------------------------------------------------------------------------------------- +!> @brief integrate state with 1st order explicit Euler method +!-------------------------------------------------------------------------------------------------- +subroutine integrateStateEuler(g,i,e) + + integer, intent(in) :: & + e, & !< element index in element loop + i, & !< integration point index in ip loop + g !< grain index in grain loop + integer :: & + p, & + c, & + sizeDotState + logical :: & + broken + + p = material_phaseAt(g,e) + c = material_phaseMemberAt(g,i,e) + + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partitionedF0, & + constitutive_mech_Fi(p)%data(1:3,1:3,c), & + crystallite_partitionedFp0, & + crystallite_subdt(g,i,e), g,i,e,p,c) + 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) + + 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(broken) return + + broken = integrateStress(g,i,e) + crystallite_converged(g,i,e) = .not. broken + +end subroutine integrateStateEuler + + +!-------------------------------------------------------------------------------------------------- +!> @brief integrate stress, state with 1st order Euler method with adaptive step size +!-------------------------------------------------------------------------------------------------- +subroutine integrateStateAdaptiveEuler(g,i,e) + + integer, intent(in) :: & + e, & !< element index in element loop + i, & !< integration point index in ip loop + g !< grain index in grain loop + integer :: & + p, & + c, & + sizeDotState + logical :: & + broken + + real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: residuum_plastic + + + p = material_phaseAt(g,e) + c = material_phaseMemberAt(g,i,e) + + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partitionedF0, & + constitutive_mech_Fi(p)%data(1:3,1:3,c), & + crystallite_partitionedFp0, & + crystallite_subdt(g,i,e), g,i,e,p,c) + if(broken) return + + sizeDotState = plasticState(p)%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) + + 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(broken) return + + broken = integrateStress(g,i,e) + if(broken) return + + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partitionedF0, & + constitutive_mech_Fi(p)%data(1:3,1:3,c), & + crystallite_partitionedFp0, & + crystallite_subdt(g,i,e), g,i,e,p,c) + 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)) + +end subroutine integrateStateAdaptiveEuler + + +!--------------------------------------------------------------------------------------------------- +!> @brief Integrate state (including stress integration) with the classic Runge Kutta method +!--------------------------------------------------------------------------------------------------- +subroutine integrateStateRK4(g,i,e) + + integer, intent(in) :: g,i,e + + real(pReal), dimension(3,3), parameter :: & + A = reshape([& + 0.5_pReal, 0.0_pReal, 0.0_pReal, & + 0.0_pReal, 0.5_pReal, 0.0_pReal, & + 0.0_pReal, 0.0_pReal, 1.0_pReal],& + shape(A)) + real(pReal), dimension(3), parameter :: & + C = [0.5_pReal, 0.5_pReal, 1.0_pReal] + 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) + +end subroutine integrateStateRK4 + + +!--------------------------------------------------------------------------------------------------- +!> @brief Integrate state (including stress integration) with the Cash-Carp method +!--------------------------------------------------------------------------------------------------- +subroutine integrateStateRKCK45(g,i,e) + + integer, intent(in) :: g,i,e + + real(pReal), dimension(5,5), parameter :: & + A = reshape([& + 1._pReal/5._pReal, .0_pReal, .0_pReal, .0_pReal, .0_pReal, & + 3._pReal/40._pReal, 9._pReal/40._pReal, .0_pReal, .0_pReal, .0_pReal, & + 3_pReal/10._pReal, -9._pReal/10._pReal, 6._pReal/5._pReal, .0_pReal, .0_pReal, & + -11._pReal/54._pReal, 5._pReal/2._pReal, -70.0_pReal/27.0_pReal, 35.0_pReal/27.0_pReal, .0_pReal, & + 1631._pReal/55296._pReal,175._pReal/512._pReal,575._pReal/13824._pReal,44275._pReal/110592._pReal,253._pReal/4096._pReal],& + shape(A)) + real(pReal), dimension(5), parameter :: & + C = [0.2_pReal, 0.3_pReal, 0.6_pReal, 1.0_pReal, 0.875_pReal] + real(pReal), dimension(6), parameter :: & + B = & + [37.0_pReal/378.0_pReal, .0_pReal, 250.0_pReal/621.0_pReal, & + 125.0_pReal/594.0_pReal, .0_pReal, 512.0_pReal/1771.0_pReal], & + DB = B - & + [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) + +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) + + + real(pReal), dimension(:,:), intent(in) :: A + real(pReal), dimension(:), intent(in) :: B, CC + 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 + integer :: & + stage, & ! stage index in integration stage loop + n, & + p, & + c, & + sizeDotState + logical :: & + broken + real(pReal), dimension(constitutive_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState + + p = material_phaseAt(g,e) + c = material_phaseMemberAt(g,i,e) + + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partitionedF0, & + constitutive_mech_Fi(p)%data(1:3,1:3,c), & + crystallite_partitionedFp0, & + crystallite_subdt(g,i,e), g,i,e,p,c) + 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) + + do n = 2, stage + sizeDotState = plasticState(p)%sizeDotState + plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) & + + 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) + + broken = integrateStress(g,i,e,CC(stage)) + if(broken) exit + + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partitionedF0, & + constitutive_mech_Fi(p)%data(1:3,1:3,c), & + crystallite_partitionedFp0, & + crystallite_subdt(g,i,e)*CC(stage), g,i,e,p,c) + if(broken) exit + + enddo + if(broken) return + + sizeDotState = plasticState(p)%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) + 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)) + + 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) + if(broken) return + + broken = integrateStress(g,i,e) + crystallite_converged(g,i,e) = .not. broken + + +end subroutine integrateStateRK + + end submodule constitutive_mech