diff --git a/code/CPFEM.f90 b/code/CPFEM.f90 index 7f375da4f..457dd12e8 100644 --- a/code/CPFEM.f90 +++ b/code/CPFEM.f90 @@ -118,6 +118,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt crystallite_Fp, & crystallite_Lp0, & crystallite_Lp, & + crystallite_dPdF0, & + crystallite_dPdF, & crystallite_Tstar0_v, & crystallite_Tstar_v use homogenization, only: homogenization_init, & @@ -242,6 +244,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt crystallite_F0 = crystallite_partionedF ! crystallite deformation (_subF is perturbed...) crystallite_Fp0 = crystallite_Fp ! crystallite plastic deformation crystallite_Lp0 = crystallite_Lp ! crystallite plastic velocity + crystallite_dPdF0 = crystallite_dPdF ! crystallite stiffness crystallite_Tstar0_v = crystallite_Tstar_v ! crystallite 2nd Piola Kirchhoff stress forall ( i = 1:homogenization_maxNgrains, & j = 1:mesh_maxNips, & diff --git a/code/IO.f90 b/code/IO.f90 index ac4ea4ef2..c818e89ac 100644 --- a/code/IO.f90 +++ b/code/IO.f90 @@ -1189,6 +1189,8 @@ endfunction case (294) msg = 'Non-positive tolerance for deformation gradient' + case (298) + msg = 'Chosen integration method does not exist' case (299) msg = 'Chosen perturbation method does not exist' diff --git a/code/constitutive.f90 b/code/constitutive.f90 index 3689916e1..6a5204c76 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -19,10 +19,12 @@ MODULE constitutive constitutive_state, & ! pointer array to current microstructure (end of converged time step) constitutive_state_backup, & ! pointer array to backed up microstructure (end of converged time step) constitutive_dotState, & ! pointer array to evolution of current microstructure - constitutive_dotState_backup, & ! pointer array to backed up evolution of current microstructure constitutive_previousDotState,& ! pointer array to previous evolution of current microstructure constitutive_previousDotState2,&! pointer array to 2nd previous evolution of current microstructure + constitutive_dotState_backup, & ! pointer array to backed up evolution of current microstructure + constitutive_RK4dotState, & ! pointer array to evolution of microstructure defined by classical Runge-Kutta method constitutive_relevantState ! relevant state values + type(p_vec), dimension(:,:,:,:), allocatable :: constitutive_RKCK45dotState ! pointer array to evolution of microstructure used by Cash-Karp Runge-Kutta method integer(pInt), dimension(:,:,:), allocatable :: constitutive_sizeDotState, & ! size of dotState array constitutive_sizeState, & ! size of state array per grain constitutive_sizePostResults ! size of postResults array per grain @@ -48,6 +50,7 @@ subroutine constitutive_init() !************************************** use prec, only: pReal,pInt use debug, only: debugger, selectiveDebugger, debug_e, debug_i, debug_g + use numerics, only: integrator, integratorStiffness use IO, only: IO_error, IO_open_file, IO_open_jobFile use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips use material @@ -58,7 +61,7 @@ subroutine constitutive_init() use constitutive_nonlocal integer(pInt), parameter :: fileunit = 200 - integer(pInt) e,i,g,p,myInstance,myNgrains + integer(pInt) e,i,g,p,s,myInstance,myNgrains integer(pInt), dimension(:,:), pointer :: thisSize character(len=64), dimension(:,:), pointer :: thisOutput logical knownConstitution @@ -71,7 +74,6 @@ subroutine constitutive_init() call constitutive_dislotwin_init(fileunit) call constitutive_nonlocal_init(fileunit) - close(fileunit) ! write description file for constitutive phase output @@ -123,12 +125,18 @@ subroutine constitutive_init() allocate(constitutive_state_backup(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) allocate(constitutive_dotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) allocate(constitutive_dotState_backup(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) - allocate(constitutive_previousDotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) - allocate(constitutive_previousDotState2(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) allocate(constitutive_relevantState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) allocate(constitutive_sizeDotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_sizeDotState = 0_pInt allocate(constitutive_sizeState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_sizeState = 0_pInt allocate(constitutive_sizePostResults(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)); constitutive_sizePostResults = 0_pInt + if (integrator == 1 .or. integratorStiffness == 1) then + allocate(constitutive_previousDotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) + allocate(constitutive_previousDotState2(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) + endif + if (integrator == 3 .or. integratorStiffness == 3) & + allocate(constitutive_RK4dotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) + if (integrator == 4 .or. integratorStiffness == 4) & + allocate(constitutive_RKCK45dotState(6,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) do e = 1,mesh_NcpElems ! loop over elements myNgrains = homogenization_Ngrains(mesh_element(3,e)) @@ -147,8 +155,17 @@ subroutine constitutive_init() allocate(constitutive_relevantState(g,i,e)%p(constitutive_j2_sizeState(myInstance))) allocate(constitutive_dotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) - allocate(constitutive_previousDotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) - allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) + if (integrator == 1 .or. integratorStiffness == 1) then + allocate(constitutive_previousDotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) + allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) + endif + if (integrator == 3 .or. integratorStiffness == 3) & + allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) + if (integrator == 4 .or. integratorStiffness == 4) then + do s = 1,6 + allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) + enddo + endif constitutive_state0(g,i,e)%p = constitutive_j2_stateInit(myInstance) constitutive_relevantState(g,i,e)%p = constitutive_j2_relevantState(myInstance) constitutive_sizeState(g,i,e) = constitutive_j2_sizeState(myInstance) @@ -164,8 +181,17 @@ subroutine constitutive_init() allocate(constitutive_relevantState(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) allocate(constitutive_dotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) - allocate(constitutive_previousDotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) - allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) + if (integrator == 1 .or. integratorStiffness == 1) then + allocate(constitutive_previousDotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) + allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) + endif + if (integrator == 3 .or. integratorStiffness == 3) & + allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) + if (integrator == 4 .or. integratorStiffness == 4) then + do s = 1,6 + allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) + enddo + endif constitutive_state0(g,i,e)%p = constitutive_phenopowerlaw_stateInit(myInstance) constitutive_relevantState(g,i,e)%p = constitutive_phenopowerlaw_relevantState(myInstance) constitutive_sizeState(g,i,e) = constitutive_phenopowerlaw_sizeState(myInstance) @@ -181,8 +207,17 @@ subroutine constitutive_init() allocate(constitutive_relevantState(g,i,e)%p(constitutive_titanmod_sizeState(myInstance))) allocate(constitutive_dotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) - allocate(constitutive_previousDotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) - allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) + if (integrator == 1 .or. integratorStiffness == 1) then + allocate(constitutive_previousDotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) + allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) + endif + if (integrator == 3 .or. integratorStiffness == 3) & + allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) + if (integrator == 4 .or. integratorStiffness == 4) then + do s = 1,6 + allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) + enddo + endif constitutive_state0(g,i,e)%p = constitutive_titanmod_stateInit(myInstance) constitutive_relevantState(g,i,e)%p = constitutive_titanmod_relevantState(myInstance) constitutive_sizeState(g,i,e) = constitutive_titanmod_sizeState(myInstance) @@ -198,8 +233,17 @@ subroutine constitutive_init() allocate(constitutive_relevantState(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance))) allocate(constitutive_dotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) - allocate(constitutive_previousDotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) - allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) + if (integrator == 1 .or. integratorStiffness == 1) then + allocate(constitutive_previousDotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) + allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) + endif + if (integrator == 3 .or. integratorStiffness == 3) & + allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) + if (integrator == 4 .or. integratorStiffness == 4) then + do s = 1,6 + allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) + enddo + endif constitutive_state0(g,i,e)%p = constitutive_dislotwin_stateInit(myInstance) constitutive_relevantState(g,i,e)%p = constitutive_dislotwin_relevantState(myInstance) constitutive_sizeState(g,i,e) = constitutive_dislotwin_sizeState(myInstance) @@ -215,8 +259,17 @@ subroutine constitutive_init() allocate(constitutive_relevantState(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance))) allocate(constitutive_dotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) - allocate(constitutive_previousDotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) - allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) + if (integrator == 1 .or. integratorStiffness == 1) then + allocate(constitutive_previousDotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) + allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) + endif + if (integrator == 3 .or. integratorStiffness == 3) & + allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) + if (integrator == 4 .or. integratorStiffness == 4) then + do s = 1,6 + allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) + enddo + endif constitutive_state0(g,i,e)%p = constitutive_nonlocal_stateInit(myInstance) constitutive_relevantState(g,i,e)%p = constitutive_nonlocal_relevantState(myInstance) constitutive_sizeState(g,i,e) = constitutive_nonlocal_sizeState(myInstance) @@ -231,7 +284,7 @@ subroutine constitutive_init() enddo enddo enddo - + constitutive_maxSizeState = maxval(constitutive_sizeState) constitutive_maxSizeDotState = maxval(constitutive_sizeDotState) constitutive_maxSizePostResults = maxval(constitutive_sizePostResults) @@ -247,7 +300,6 @@ subroutine constitutive_init() write(6,'(a32,x,7(i5,x))') 'constitutive_state: ', shape(constitutive_state) write(6,'(a32,x,7(i5,x))') 'constitutive_relevantState: ', shape(constitutive_relevantState) write(6,'(a32,x,7(i5,x))') 'constitutive_dotState: ', shape(constitutive_dotState) - write(6,'(a32,x,7(i5,x))') 'constitutive_previousDotState:', shape(constitutive_previousDotState) write(6,'(a32,x,7(i5,x))') 'constitutive_sizeState: ', shape(constitutive_sizeState) write(6,'(a32,x,7(i5,x))') 'constitutive_sizeDotState: ', shape(constitutive_sizeDotState) write(6,'(a32,x,7(i5,x))') 'constitutive_sizePostResults: ', shape(constitutive_sizePostResults) @@ -448,7 +500,8 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, Temperature, ip call constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el) case (constitutive_nonlocal_label) - call constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, Temperature, constitutive_state, ipc, ip, el) + call constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, Temperature, constitutive_state, & + constitutive_relevantState, ipc, ip, el) end select @@ -478,11 +531,16 @@ use mesh, only: mesh_NcpElems, & use material, only: phase_constitution, & material_phase, & homogenization_maxNgrains -use constitutive_j2 -use constitutive_phenopowerlaw -use constitutive_titanmod -use constitutive_dislotwin -use constitutive_nonlocal +use constitutive_j2, only: constitutive_j2_dotState, & + constitutive_j2_label +use constitutive_phenopowerlaw, only: constitutive_phenopowerlaw_dotState, & + constitutive_phenopowerlaw_label +use constitutive_titanmod, only: constitutive_titanmod_dotState, & + constitutive_titanmod_label +use constitutive_dislotwin, only: constitutive_dislotwin_dotState, & + constitutive_dislotwin_label +use constitutive_nonlocal, only: constitutive_nonlocal_dotState, & + constitutive_nonlocal_label implicit none @@ -548,40 +606,59 @@ function constitutive_dotTemperature(Tstar_v,Temperature,ipc,ip,el) !* OUTPUT: * !* - constitutive_dotTemperature : evolution of temperature * !********************************************************************* - use prec, only: pReal,pInt - use material, only: phase_constitution,material_phase - use constitutive_j2 - use constitutive_phenopowerlaw - use constitutive_titanmod - use constitutive_dislotwin - use constitutive_nonlocal - implicit none +use prec, only: pReal,pInt +use debug, only: debug_cumDotTemperatureCalls, & + debug_cumDotTemperatureTicks +use material, only: phase_constitution,material_phase +use constitutive_j2, only: constitutive_j2_dotTemperature, & + constitutive_j2_label +use constitutive_phenopowerlaw, only: constitutive_phenopowerlaw_dotTemperature, & + constitutive_phenopowerlaw_label +use constitutive_titanmod, only: constitutive_titanmod_dotTemperature, & + constitutive_titanmod_label +use constitutive_dislotwin, only: constitutive_dislotwin_dotTemperature, & + constitutive_dislotwin_label +use constitutive_nonlocal, only: constitutive_nonlocal_dotTemperature, & + constitutive_nonlocal_label +implicit none !* Definition of variables - integer(pInt) ipc,ip,el - real(pReal) Temperature - real(pReal) constitutive_dotTemperature - real(pReal), dimension(6) :: Tstar_v +integer(pInt) ipc, ip, el +real(pReal) Temperature +real(pReal) constitutive_dotTemperature +real(pReal), dimension(6) :: Tstar_v +integer(pLongInt) tick, & + tock, & + tickrate, & + maxticks - select case (phase_constitution(material_phase(ipc,ip,el))) - - case (constitutive_j2_label) - constitutive_dotTemperature = constitutive_j2_dotTemperature(Tstar_v,Temperature,constitutive_state,ipc,ip,el) - - case (constitutive_phenopowerlaw_label) - constitutive_dotTemperature = constitutive_phenopowerlaw_dotTemperature(Tstar_v,Temperature,constitutive_state,ipc,ip,el) +call system_clock(count=tick,count_rate=tickrate,count_max=maxticks) - case (constitutive_titanmod_label) - constitutive_dotTemperature = constitutive_titanmod_dotTemperature(Tstar_v,Temperature,constitutive_state,ipc,ip,el) - - case (constitutive_dislotwin_label) - constitutive_dotTemperature = constitutive_dislotwin_dotTemperature(Tstar_v,Temperature,constitutive_state,ipc,ip,el) - - case (constitutive_nonlocal_label) - constitutive_dotTemperature = constitutive_nonlocal_dotTemperature(Tstar_v,Temperature,constitutive_state,ipc,ip,el) - - end select - return +select case (phase_constitution(material_phase(ipc,ip,el))) + + case (constitutive_j2_label) + constitutive_dotTemperature = constitutive_j2_dotTemperature(Tstar_v,Temperature,constitutive_state,ipc,ip,el) + + case (constitutive_phenopowerlaw_label) + constitutive_dotTemperature = constitutive_phenopowerlaw_dotTemperature(Tstar_v,Temperature,constitutive_state,ipc,ip,el) + + case (constitutive_titanmod_label) + constitutive_dotTemperature = constitutive_titanmod_dotTemperature(Tstar_v,Temperature,constitutive_state,ipc,ip,el) + + case (constitutive_dislotwin_label) + constitutive_dotTemperature = constitutive_dislotwin_dotTemperature(Tstar_v,Temperature,constitutive_state,ipc,ip,el) + + case (constitutive_nonlocal_label) + constitutive_dotTemperature = constitutive_nonlocal_dotTemperature(Tstar_v,Temperature,constitutive_state,ipc,ip,el) + +end select + +call system_clock(count=tock,count_rate=tickrate,count_max=maxticks) +debug_cumDotTemperatureCalls = debug_cumDotTemperatureCalls + 1_pInt +debug_cumDotTemperatureTicks = debug_cumDotTemperatureTicks + tock-tick +if (tock < tick) debug_cumDotTemperatureTicks = debug_cumDotTemperatureTicks + maxticks + +return endfunction diff --git a/code/constitutive_nonlocal.f90 b/code/constitutive_nonlocal.f90 index 9104eb6ea..c5c485990 100644 --- a/code/constitutive_nonlocal.f90 +++ b/code/constitutive_nonlocal.f90 @@ -257,7 +257,6 @@ constitutive_nonlocal_interactionSlipSlip = 0.0_pReal constitutive_nonlocal_dLowerEdgePerSlipFamily = 0.0_pReal constitutive_nonlocal_dLowerScrewPerSlipFamily = 0.0_pReal - !*** readout data from material.config file rewind(file) @@ -973,7 +972,7 @@ endsubroutine !********************************************************************* !* calculates kinetics * !********************************************************************* -subroutine constitutive_nonlocal_kinetics(Tstar_v, Temperature, state, g, ip, el) +subroutine constitutive_nonlocal_kinetics(Tstar_v, Temperature, state, g, ip, el, dv_dtau) use prec, only: pReal, & pInt, & @@ -1003,6 +1002,8 @@ type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), in real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola-Kirchhoff stress in Mandel notation !*** output variables +real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))), & + intent(out), optional :: dv_dtau ! velocity derivative with respect to resolved shear stress !*** local variables integer(pInt) myInstance, & ! current instance of this constitution @@ -1013,31 +1014,42 @@ integer(pInt) myInstance, & ! curren real(pReal), dimension(6) :: Tdislocation_v ! dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: & tauThreshold, & ! threshold shear stress - tau ! resolved shear stress + tau, & ! resolved shear stress + rhoForest ! forest dislocation density +real(pReal) boltzmannProbability myInstance = phase_constitutionInstance(material_phase(g,ip,el)) myStructure = constitutive_nonlocal_structure(myInstance) ns = constitutive_nonlocal_totalNslip(myInstance) +rhoForest = state(g,ip,el)%p(10*ns+1:11*ns) tauThreshold = state(g,ip,el)%p(11*ns+1:12*ns) Tdislocation_v = state(g,ip,el)%p(12*ns+1:12*ns+6) tau = 0.0_pReal constitutive_nonlocal_v(:,:,g,ip,el) = 0.0_pReal +if (present(dv_dtau)) dv_dtau = 0.0_pReal -do s = 1,ns - if ((tauThreshold(s) > 0.0_pReal) .and. (Temperature > 0.0_pReal)) then +if ( Temperature > 0.0_pReal ) then + do s = 1,ns tau(s) = math_mul6x6( Tstar_v + Tdislocation_v, & - lattice_Sslip_v(:,constitutive_nonlocal_slipSystemLattice(s,myInstance),myStructure) ) - - constitutive_nonlocal_v(s,:,g,ip,el) = constitutive_nonlocal_v0PerSlipSystem(s,myInstance) & - * exp( - constitutive_nonlocal_Q0(myInstance) / ( kB * Temperature) * (1.0_pReal - (abs(tau(s))/tauThreshold(s)) ) ) & - * sign(1.0_pReal,tau(s)) - - endif -enddo + lattice_Sslip_v(:,constitutive_nonlocal_slipSystemLattice(s,myInstance),myStructure) ) + + if ( abs(tau(s)) > 0.0_pReal ) then + boltzmannProbability = dexp( -constitutive_nonlocal_Q0(myInstance) * dsqrt(rhoForest(s)) / ( abs(tau(s)) * kB * Temperature) ) + constitutive_nonlocal_v(s,:,g,ip,el) = constitutive_nonlocal_v0PerSlipSystem(s,myInstance) & + / ( constitutive_nonlocal_burgersPerSlipSystem(s,myInstance) * dsqrt(rhoForest(s)) ) & + * boltzmannProbability * sign(1.0_pReal,tau(s)) + + if (present(dv_dtau)) & + dv_dtau(s) = constitutive_nonlocal_Q0(myInstance) * constitutive_nonlocal_v0PerSlipSystem(s,myInstance) & + / ( constitutive_nonlocal_burgersPerSlipSystem(s,myInstance) * tau(s)**2.0_pReal * kB * Temperature ) & + * boltzmannProbability + endif + enddo +endif if (verboseDebugger .and. selectiveDebugger) then !$OMP CRITICAL (write2out) @@ -1058,7 +1070,7 @@ endsubroutine !********************************************************************* !* calculates plastic velocity gradient and its tangent * !********************************************************************* -subroutine constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar99, Tstar_v, Temperature, state, g, ip, el) +subroutine constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar99, Tstar_v, Temperature, state, relevantState, g, ip, el) use prec, only: pReal, & pInt, & @@ -1085,7 +1097,8 @@ integer(pInt), intent(in) :: g, & ! curren el ! current element number real(pReal), intent(in) :: Temperature ! temperature type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & - state ! microstructural state + state, & ! microstructural state + relevantState ! relevant microstructural state real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola-Kirchhoff stress in Mandel notation !*** output variables @@ -1111,7 +1124,9 @@ real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstan real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: & tauThreshold, & ! threshold shear stress gdotTotal, & ! shear rate - dgdotTotal_dtau ! derivative of the shear rate with respect to the shear stress + dv_dtau, & ! velocity derivative with respect to the shear stress + dgdotTotal_dtau, & ! derivative of the shear rate with respect to the shear stress + rhoForest ! forest dislocation density !*** initialize local variables @@ -1131,17 +1146,18 @@ forall (t = 1:8) & forall (s = 1:ns, t = 5:8, rhoSgl(s,t) * constitutive_nonlocal_v(s,t-4,g,ip,el) < 0.0_pReal) & ! contribution of used rho for changing sign of v rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t)) +rhoForest = state(g,ip,el)%p(10*ns+1:11*ns) tauThreshold = state(g,ip,el)%p(11*ns+1:12*ns) -call constitutive_nonlocal_kinetics(Tstar_v, Temperature, state, g, ip, el) ! update dislocation velocity +call constitutive_nonlocal_kinetics(Tstar_v, Temperature, state, g, ip, el, dv_dtau) ! update dislocation velocity !*** Calculation of gdot and its tangent -forall (t = 1:4 ) & - gdot(:,t) = rhoSgl(:,t) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) * constitutive_nonlocal_v(:,t,g,ip,el) +forall (s = 1:ns, t = 1:4, rhoSgl(s,t) > relevantState(g,ip,el)%p((t-1)*ns+s)) & ! no shear rate contribution for densities below relevant state + gdot(s,t) = rhoSgl(s,t) * constitutive_nonlocal_burgersPerSlipSystem(s,myInstance) * constitutive_nonlocal_v(s,t,g,ip,el) gdotTotal = sum(gdot,2) -dgdotTotal_dtau = abs(gdotTotal) * constitutive_nonlocal_Q0(myInstance) / ( kB * Temperature * tauThreshold ) +dgdotTotal_dtau = sum(rhoSgl,2) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) * dv_dtau !*** Calculation of Lp and its tangent @@ -1216,7 +1232,8 @@ use lattice, only: lattice_Sslip, & lattice_sn, & lattice_st, & lattice_maxNslipFamily, & - lattice_NslipSystem + lattice_NslipSystem +use FEsolving, only:theInc implicit none !*** input variables @@ -1352,16 +1369,16 @@ if (timestep <= 0.0_pReal) then return endif -if (any(constitutive_nonlocal_v(:,:,g,ip,el)*timestep > mesh_ipVolume(ip,el)**(1.0_pReal/3.0_pReal))) then ! if timestep is too large,... - dotState(1,ip,el)%p(1:10*ns) = NaN ! ...assign NaN and enforce a cutback - if (verboseDebugger) then - !$OMP CRITICAL (write2out) - write(6,*) 'exceeded maximum allowed dislocation velocity at ',g,ip,el - write(6,*) - !$OMPEND CRITICAL (write2out) - endif - return -endif +!if (any(constitutive_nonlocal_v(:,:,g,ip,el)*timestep > mesh_ipVolume(ip,el)**(1.0_pReal/3.0_pReal))) then ! if timestep is too large,... +! dotState(1,ip,el)%p(1:10*ns) = NaN ! ...assign NaN and enforce a cutback +! if (verboseDebugger) then +! !$OMP CRITICAL (write2out) +! write(6,*) 'exceeded maximum allowed dislocation velocity at ',g,ip,el +! write(6,*) +! !$OMPEND CRITICAL (write2out) +! endif +! return +!endif !**************************************************************************** @@ -1603,20 +1620,23 @@ forall (t = 1:10) & + rhoDotAthermalAnnihilation(:,t) & + rhoDotRemobilization(:,t) & + rhoDotMultiplication(:,t) & - + rhoDotThermalAnnihilation(:,t) -! + rhoDotDipole2SingleStressChange(:,t) & + + rhoDotThermalAnnihilation(:,t) +! + rhoDotDipole2SingleStressChange(:,t) ! + rhoDotSingle2DipoleStressChange(:,t) -dotState(g,ip,el)%p(1:10*ns) = dotState(g,ip,el)%p(1:10*ns) + reshape(rhoDot,(/10*ns/)) +dotState(g,ip,el)%p(1:10*ns) = reshape(rhoDot,(/10*ns/)) -do i = 1,4*ns +do i = 1,10*ns + if (i > 4*ns .and. i <= 8*ns) & ! skip immobile densities + continue if (previousState(g,ip,el)%p(i) + dotState(g,ip,el)%p(i)*timestep < 0.0_pReal) then ! if single mobile densities become negative... if (previousState(g,ip,el)%p(i) < relevantState(g,ip,el)%p(i)) then ! ... and density is already below relevance... - dotState(g,ip,el)%p(i) = 0.0_pReal ! ... set dotState to zero + dotState(g,ip,el)%p(i) = - previousState(g,ip,el)%p(i) / timestep ! ... set dotState to zero else ! ... otherwise... if (verboseDebugger) then !$OMP CRITICAL (write2out) - write(6,*) 'negative dislocation density at ',g,ip,el + write(6,*) 'negative dislocation density: enforced cutback at ',g,ip,el,i + write(6,'(a12,x,e15.8)') 'dotState was',dotState(g,ip,el)%p(i) write(6,*) !$OMPEND CRITICAL (write2out) endif @@ -1671,7 +1691,7 @@ disorientationAxisAngle = math_QuaternionToAxisAngle(disorientation) disorientationAxis = disorientationAxisAngle(1:3) disorientationAngle = disorientationAxisAngle(4) -if (disorientationAngle < 3.0_pReal) then +if (disorientationAngle < 5.0_pReal) then constitutive_nonlocal_transmissivity = 1.0_pReal else constitutive_nonlocal_transmissivity = 0.5_pReal diff --git a/code/crystallite.f90 b/code/crystallite.f90 index a729363d5..dcf308b82 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -10,7 +10,7 @@ !* - _stressAndItsTangent * !* - _postResults * !*************************************** - + MODULE crystallite use prec, only: pReal, pInt @@ -32,10 +32,11 @@ real(pReal), dimension (:,:,:), allocatable :: & crystallite_subdt, & ! substepped time increment of each grain crystallite_subFrac, & ! already calculated fraction of increment crystallite_subStep, & ! size of next integration step + crystallite_statedamper, & ! damping for state update crystallite_Temperature, & ! Temp of each grain crystallite_partionedTemperature0, & ! Temp of each grain at start of homog inc crystallite_subTemperature0, & ! Temp of each grain at start of crystallite inc - crystallite_statedamper ! damping for state update + crystallite_dotTemperature ! evolution of Temperature of each grain real(pReal), dimension (:,:,:,:), allocatable :: & crystallite_Tstar_v, & ! current 2nd Piola-Kirchhoff stress vector (end of converged time step) crystallite_Tstar0_v, & ! 2nd Piola-Kirchhoff stress vector at start of FE inc @@ -63,12 +64,14 @@ real(pReal), dimension (:,:,:,:,:), allocatable :: & crystallite_P, & ! 1st Piola-Kirchhoff stress per grain crystallite_disorientation ! disorientation between two neighboring ips (only calculated for single grain IPs) real(pReal), dimension (:,:,:,:,:,:,:), allocatable :: & - crystallite_dPdF, & ! individual dPdF per grain + crystallite_dPdF, & ! current individual dPdF per grain (end of converged time step) + crystallite_dPdF0, & ! individual dPdF per grain at start of FE inc + crystallite_partioneddPdF0, & ! individual dPdF per grain at start of homog inc crystallite_fallbackdPdF ! dPdF fallback for non-converged grains (elastic prediction) logical, dimension (:,:,:), allocatable :: & crystallite_localConstitution, & ! indicates this grain to have purely local constitutive law crystallite_requested, & ! flag to request crystallite calculation - crystallite_todo, & ! flag to indicate ongoing calculation + crystallite_todo, & ! flag to indicate need for further computation crystallite_converged, & ! convergence flag crystallite_stateConverged, & ! flag indicating convergence of state crystallite_temperatureConverged ! flag indicating convergence of temperature @@ -85,6 +88,10 @@ subroutine crystallite_init(Temperature) pReal use debug, only: debug_info, & debug_reset + use numerics, only: integrator, & + integratorStiffness, & + subStepSizeCryst, & + stepIncreaseCryst use math, only: math_I3, & math_EulerToR use FEsolving, only: FEsolving_execElem, & @@ -95,7 +102,7 @@ subroutine crystallite_init(Temperature) mesh_maxNipNeighbors use IO use material - use lattice, only: lattice_symmetryType + use lattice, only: lattice_symmetryTypes use constitutive_phenopowerlaw, only: constitutive_phenopowerlaw_label, & constitutive_phenopowerlaw_structure use constitutive_titanmod, only: constitutive_titanmod_label, & @@ -143,47 +150,50 @@ subroutine crystallite_init(Temperature) nMax = mesh_maxNipNeighbors allocate(crystallite_Temperature(gMax,iMax,eMax)); crystallite_Temperature = Temperature + allocate(crystallite_partionedTemperature0(gMax,iMax,eMax)); crystallite_partionedTemperature0 = 0.0_pReal + allocate(crystallite_subTemperature0(gMax,iMax,eMax)); crystallite_subTemperature0 = 0.0_pReal + allocate(crystallite_dotTemperature(gMax,iMax,eMax)); crystallite_dotTemperature = 0.0_pReal + allocate(crystallite_Tstar0_v(6,gMax,iMax,eMax)); crystallite_Tstar0_v = 0.0_pReal + allocate(crystallite_partionedTstar0_v(6,gMax,iMax,eMax)); crystallite_partionedTstar0_v = 0.0_pReal + allocate(crystallite_subTstar0_v(6,gMax,iMax,eMax)); crystallite_subTstar0_v = 0.0_pReal + allocate(crystallite_Tstar_v(6,gMax,iMax,eMax)); crystallite_Tstar_v = 0.0_pReal allocate(crystallite_P(3,3,gMax,iMax,eMax)); crystallite_P = 0.0_pReal - allocate(crystallite_Fe(3,3,gMax,iMax,eMax)); crystallite_Fe = 0.0_pReal + allocate(crystallite_F0(3,3,gMax,iMax,eMax)); crystallite_F0 = 0.0_pReal + allocate(crystallite_partionedF0(3,3,gMax,iMax,eMax)); crystallite_partionedF0 = 0.0_pReal + allocate(crystallite_partionedF(3,3,gMax,iMax,eMax)); crystallite_partionedF = 0.0_pReal + allocate(crystallite_subF0(3,3,gMax,iMax,eMax)); crystallite_subF0 = 0.0_pReal + allocate(crystallite_subF(3,3,gMax,iMax,eMax)); crystallite_subF = 0.0_pReal + allocate(crystallite_Fp0(3,3,gMax,iMax,eMax)); crystallite_Fp0 = 0.0_pReal + allocate(crystallite_partionedFp0(3,3,gMax,iMax,eMax)); crystallite_partionedFp0 = 0.0_pReal + allocate(crystallite_subFp0(3,3,gMax,iMax,eMax)); crystallite_subFp0 = 0.0_pReal allocate(crystallite_Fp(3,3,gMax,iMax,eMax)); crystallite_Fp = 0.0_pReal allocate(crystallite_invFp(3,3,gMax,iMax,eMax)); crystallite_invFp = 0.0_pReal - allocate(crystallite_Lp(3,3,gMax,iMax,eMax)); crystallite_Lp = 0.0_pReal - allocate(crystallite_Tstar_v(6,gMax,iMax,eMax)); crystallite_Tstar_v = 0.0_pReal - allocate(crystallite_F0(3,3,gMax,iMax,eMax)); crystallite_F0 = 0.0_pReal - allocate(crystallite_Fp0(3,3,gMax,iMax,eMax)); crystallite_Fp0 = 0.0_pReal + allocate(crystallite_Fe(3,3,gMax,iMax,eMax)); crystallite_Fe = 0.0_pReal allocate(crystallite_Lp0(3,3,gMax,iMax,eMax)); crystallite_Lp0 = 0.0_pReal - allocate(crystallite_Tstar0_v(6,gMax,iMax,eMax)); crystallite_Tstar0_v = 0.0_pReal - allocate(crystallite_partionedTemperature0(gMax,iMax,eMax)); crystallite_partionedTemperature0 = 0.0_pReal - allocate(crystallite_partionedF(3,3,gMax,iMax,eMax)); crystallite_partionedF = 0.0_pReal - allocate(crystallite_partionedF0(3,3,gMax,iMax,eMax)); crystallite_partionedF0 = 0.0_pReal - allocate(crystallite_partionedFp0(3,3,gMax,iMax,eMax)); crystallite_partionedFp0 = 0.0_pReal allocate(crystallite_partionedLp0(3,3,gMax,iMax,eMax)); crystallite_partionedLp0 = 0.0_pReal - allocate(crystallite_partionedTstar0_v(6,gMax,iMax,eMax)); crystallite_partionedTstar0_v = 0.0_pReal - allocate(crystallite_subTemperature0(gMax,iMax,eMax)); crystallite_subTemperature0 = 0.0_pReal - allocate(crystallite_symmetryID(gMax,iMax,eMax)); crystallite_symmetryID = 0.0_pReal !NEW - allocate(crystallite_subF(3,3,gMax,iMax,eMax)); crystallite_subF = 0.0_pReal - allocate(crystallite_subF0(3,3,gMax,iMax,eMax)); crystallite_subF0 = 0.0_pReal - allocate(crystallite_subFp0(3,3,gMax,iMax,eMax)); crystallite_subFp0 = 0.0_pReal allocate(crystallite_subLp0(3,3,gMax,iMax,eMax)); crystallite_subLp0 = 0.0_pReal - allocate(crystallite_orientation(4,gMax,iMax,eMax)); crystallite_orientation = 0.0_pReal - allocate(crystallite_orientation0(4,gMax,iMax,eMax)); crystallite_orientation = 0.0_pReal - allocate(crystallite_rotation(4,gMax,iMax,eMax)); crystallite_rotation = 0.0_pReal - allocate(crystallite_disorientation(4,nMax,gMax,iMax,eMax)); crystallite_disorientation = 0.0_pReal - allocate(crystallite_subTstar0_v(6,gMax,iMax,eMax)); crystallite_subTstar0_v = 0.0_pReal + allocate(crystallite_Lp(3,3,gMax,iMax,eMax)); crystallite_Lp = 0.0_pReal allocate(crystallite_dPdF(3,3,3,3,gMax,iMax,eMax)); crystallite_dPdF = 0.0_pReal + allocate(crystallite_dPdF0(3,3,3,3,gMax,iMax,eMax)); crystallite_dPdF0 = 0.0_pReal + allocate(crystallite_partioneddPdF0(3,3,3,3,gMax,iMax,eMax)); crystallite_partioneddPdF0 = 0.0_pReal allocate(crystallite_fallbackdPdF(3,3,3,3,gMax,iMax,eMax)); crystallite_fallbackdPdF = 0.0_pReal allocate(crystallite_dt(gMax,iMax,eMax)); crystallite_dt = 0.0_pReal allocate(crystallite_subdt(gMax,iMax,eMax)); crystallite_subdt = 0.0_pReal allocate(crystallite_subFrac(gMax,iMax,eMax)); crystallite_subFrac = 0.0_pReal allocate(crystallite_subStep(gMax,iMax,eMax)); crystallite_subStep = 0.0_pReal allocate(crystallite_statedamper(gMax,iMax,eMax)); crystallite_statedamper = 1.0_pReal + allocate(crystallite_symmetryID(gMax,iMax,eMax)); crystallite_symmetryID = 0.0_pReal !NEW + allocate(crystallite_orientation(4,gMax,iMax,eMax)); crystallite_orientation = 0.0_pReal + allocate(crystallite_orientation0(4,gMax,iMax,eMax)); crystallite_orientation = 0.0_pReal + allocate(crystallite_rotation(4,gMax,iMax,eMax)); crystallite_rotation = 0.0_pReal + allocate(crystallite_disorientation(4,nMax,gMax,iMax,eMax)); crystallite_disorientation = 0.0_pReal allocate(crystallite_localConstitution(gMax,iMax,eMax)); crystallite_localConstitution = .true. allocate(crystallite_requested(gMax,iMax,eMax)); crystallite_requested = .false. + allocate(crystallite_todo(gMax,iMax,eMax)); crystallite_todo = .false. allocate(crystallite_converged(gMax,iMax,eMax)); crystallite_converged = .true. allocate(crystallite_stateConverged(gMax,iMax,eMax)); crystallite_stateConverged = .false. allocate(crystallite_temperatureConverged(gMax,iMax,eMax)); crystallite_temperatureConverged = .false. - allocate(crystallite_todo(gMax,iMax,eMax)); crystallite_todo = .true. - + allocate(crystallite_output(maxval(crystallite_Noutput), & material_Ncrystallite)) ; crystallite_output = '' allocate(crystallite_sizePostResults(material_Ncrystallite)) ; crystallite_sizePostResults = 0_pInt @@ -302,15 +312,15 @@ subroutine crystallite_init(Temperature) myStructure = -1_pInt ! does this happen for j2 material? end select if (myStructure>0_pInt) then - crystallite_symmetryID(g,i,e) = lattice_symmetryType(myStructure) ! structure = 1(fcc) or 2(bcc) => 1; 3+(hex)=>2 + crystallite_symmetryID(g,i,e)=lattice_symmetryTypes(myStructure) ! structure = 1(fcc) or 2(bcc) => 1; 3(hex)=>2 endif enddo enddo enddo !$OMPEND PARALLEL DO - + call crystallite_orientations() - crystallite_orientation0=crystallite_orientation ! Store initial orientations for calculation of grain rotations + crystallite_orientation0 = crystallite_orientation ! Store initial orientations for calculation of grain rotations call crystallite_stressAndItsTangent(.true.) ! request elastic answers crystallite_fallbackdPdF = crystallite_dPdF ! use initial elastic stiffness as fallback @@ -318,6 +328,7 @@ subroutine crystallite_init(Temperature) ! *** Output to MARC output file *** !$OMP CRITICAL (write2out) write(6,'(a35,x,7(i5,x))') 'crystallite_Temperature: ', shape(crystallite_Temperature) + write(6,'(a35,x,7(i5,x))') 'crystallite_dotTemperature: ', shape(crystallite_dotTemperature) write(6,'(a35,x,7(i5,x))') 'crystallite_Fe: ', shape(crystallite_Fe) write(6,'(a35,x,7(i5,x))') 'crystallite_Fp: ', shape(crystallite_Fp) write(6,'(a35,x,7(i5,x))') 'crystallite_Lp: ', shape(crystallite_Lp) @@ -341,6 +352,8 @@ subroutine crystallite_init(Temperature) write(6,'(a35,x,7(i5,x))') 'crystallite_partionedTstar0_v: ', shape(crystallite_partionedTstar0_v) write(6,'(a35,x,7(i5,x))') 'crystallite_subTstar0_v: ', shape(crystallite_subTstar0_v) write(6,'(a35,x,7(i5,x))') 'crystallite_dPdF: ', shape(crystallite_dPdF) + write(6,'(a35,x,7(i5,x))') 'crystallite_dPdF0: ', shape(crystallite_dPdF0) + write(6,'(a35,x,7(i5,x))') 'crystallite_partioneddPdF0: ', shape(crystallite_partioneddPdF0) write(6,'(a35,x,7(i5,x))') 'crystallite_fallbackdPdF: ', shape(crystallite_fallbackdPdF) write(6,'(a35,x,7(i5,x))') 'crystallite_orientation: ', shape(crystallite_orientation) write(6,'(a35,x,7(i5,x))') 'crystallite_orientation0: ', shape(crystallite_orientation0) @@ -350,6 +363,7 @@ subroutine crystallite_init(Temperature) write(6,'(a35,x,7(i5,x))') 'crystallite_subdt: ', shape(crystallite_subdt) write(6,'(a35,x,7(i5,x))') 'crystallite_subFrac: ', shape(crystallite_subFrac) write(6,'(a35,x,7(i5,x))') 'crystallite_subStep: ', shape(crystallite_subStep) + write(6,'(a35,x,7(i5,x))') 'crystallite_stateDamper: ', shape(crystallite_stateDamper) write(6,'(a35,x,7(i5,x))') 'crystallite_localConstitution: ', shape(crystallite_localConstitution) write(6,'(a35,x,7(i5,x))') 'crystallite_requested: ', shape(crystallite_requested) write(6,'(a35,x,7(i5,x))') 'crystallite_todo: ', shape(crystallite_todo) @@ -385,26 +399,24 @@ subroutine crystallite_stressAndItsTangent(updateJaco) stepIncreaseCryst, & pert_Fg, & pert_method, & - nState, & nCryst, & - iJacoStiffness + iJacoStiffness, & + integratorStiffness, & + integrator use debug, only: debugger, & selectiveDebugger, & verboseDebugger, & debug_e, & debug_i, & debug_g, & - debug_CrystalliteLoopDistribution, & - debug_CrystalliteStateLoopDistribution, & - debug_StiffnessStateLoopDistribution + debug_CrystalliteLoopDistribution use IO, only: IO_warning use math, only: math_inv3x3, & math_mul33x33, & math_mul66x6, & math_Mandel6to33, & math_Mandel33to6, & - math_I3, & - math_Plain3333to99 + math_I3 use FEsolving, only: FEsolving_execElem, & FEsolving_execIP, & theInc, & @@ -425,8 +437,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) constitutive_homogenizedC, & constitutive_dotState, & constitutive_dotState_backup, & - constitutive_previousDotState, & - constitutive_previousDotState2, & constitutive_collectDotState, & constitutive_dotTemperature, & constitutive_microstructure @@ -440,33 +450,13 @@ subroutine crystallite_stressAndItsTangent(updateJaco) !*** local variables ***! real(pReal) myTemperature, & ! local copy of the temperature - myPert ! perturbation with correct sign + myPert, & ! perturbation with correct sign + formerSubStep real(pReal), dimension(3,3) :: invFp, & ! inverse of the plastic deformation gradient Fe_guess, & ! guess for elastic deformation gradient Tstar ! 2nd Piola-Kirchhoff stress tensor - integer(pInt) NiterationCrystallite, & ! number of iterations in crystallite loop - NiterationState ! number of iterations in state loop - integer(pInt) e, ee, & ! element index - i, ii, & ! integration point index - g, gg, & ! grain index - k, & - l, & - perturbation , & ! loop counter for forward,backward perturbation mode - comp, & - myNgrains, & - mySizeState, & - mySizeDotState - integer(pInt), dimension(2,9,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & - kl - logical onTrack, & ! flag indicating whether we are still on track - temperatureConverged, & ! flag indicating if temperature converged - stateConverged, & ! flag indicating if state converged - converged ! flag indicating if iteration converged real(pReal), dimension(9,9) :: dPdF99 real(pReal), dimension(3,3,3,3,2) :: dPdF_perturbation - real(pReal) dot_prod12, & - dot_prod22, & - formerSubStep real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & F_backup, & Fp_backup, & @@ -478,44 +468,44 @@ subroutine crystallite_stressAndItsTangent(updateJaco) Tstar_v_backup real(pReal), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & Temperature_backup + integer(pInt) NiterationCrystallite, & ! number of iterations in crystallite loop + e, & ! element index + i, & ! integration point index + g, & ! grain index + k, & + l, & + perturbation , & ! loop counter for forward,backward perturbation mode + myNgrains, & + mySizeState, & + mySizeDotState logical, dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & - ConvergenceFlag_backup - logical, dimension(3,3) :: mask + convergenceFlag_backup logical forceLocalStiffnessCalculation ! flag indicating that stiffness calculation is always done locally forceLocalStiffnessCalculation = .true. - ! ------ initialize to starting condition ------ + ! --+>> INITIALIZE TO STARTING CONDITION <<+-- -!$OMP CRITICAL (write2out) -! write (6,*) -! write (6,*) 'Crystallite request from Materialpoint' -! write (6,'(a,/,(f12.7,x))') 'crystallite_partionedTemp0 of 1 1 1' ,crystallite_partionedTemperature0(1,1,1) -! write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedF0 of 1 1 1' ,crystallite_partionedF0(1:3,:,1,1,1) -! write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedFp0 of 1 1 1' ,crystallite_partionedFp0(1:3,:,1,1,1) -! write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedF of 1 1 1' ,crystallite_partionedF(1:3,:,1,1,1) -! write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedLp0 of 1 1 1' ,crystallite_partionedLp0(1:3,:,1,1,1) -!$OMPEND CRITICAL (write2out) - crystallite_subStep = 0.0_pReal !$OMP PARALLEL DO - do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed + do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed myNgrains = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains - if (crystallite_requested(g,i,e)) then ! initialize restoration point of ... - crystallite_subTemperature0(g,i,e) = crystallite_partionedTemperature0(g,i,e) ! ...temperature - constitutive_subState0(g,i,e)%p = constitutive_partionedState0(g,i,e)%p ! ...microstructure - crystallite_subFp0(:,:,g,i,e) = crystallite_partionedFp0(:,:,g,i,e) ! ...plastic def grad - crystallite_subLp0(:,:,g,i,e) = crystallite_partionedLp0(:,:,g,i,e) ! ...plastic velocity grad - crystallite_subF0(:,:,g,i,e) = crystallite_partionedF0(:,:,g,i,e) ! ...def grad - crystallite_subTstar0_v(:,g,i,e) = crystallite_partionedTstar0_v(:,g,i,e) !...2nd PK stress + if (crystallite_requested(g,i,e)) then ! initialize restoration point of ... + crystallite_subTemperature0(g,i,e) = crystallite_partionedTemperature0(g,i,e) ! ...temperature + constitutive_subState0(g,i,e)%p = constitutive_partionedState0(g,i,e)%p ! ...microstructure + crystallite_subFp0(:,:,g,i,e) = crystallite_partionedFp0(:,:,g,i,e) ! ...plastic def grad + crystallite_subLp0(:,:,g,i,e) = crystallite_partionedLp0(:,:,g,i,e) ! ...plastic velocity grad + crystallite_dPdF0(:,:,:,:,g,i,e) = crystallite_partioneddPdF0(:,:,:,:,g,i,e) ! ...stiffness + crystallite_subF0(:,:,g,i,e) = crystallite_partionedF0(:,:,g,i,e) ! ...def grad + crystallite_subTstar0_v(:,g,i,e) = crystallite_partionedTstar0_v(:,g,i,e) !...2nd PK stress crystallite_subFrac(g,i,e) = 0.0_pReal crystallite_subStep(g,i,e) = 1.0_pReal/subStepSizeCryst crystallite_todo(g,i,e) = .true. - crystallite_converged(g,i,e) = .false. ! pretend failed step of twice the required size + crystallite_converged(g,i,e) = .false. ! pretend failed step of twice the required size endif enddo enddo @@ -523,18 +513,20 @@ subroutine crystallite_stressAndItsTangent(updateJaco) !$OMPEND PARALLEL DO - ! --+>> crystallite loop <<+-- + ! --+>> CRYSTALLITE CUTBACK LOOP <<+-- NiterationCrystallite = 0_pInt - - do while (any(crystallite_subStep(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) > subStepMinCryst)) ! cutback loop for crystallites + do while (any(crystallite_subStep(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) > subStepMinCryst)) ! cutback loop for crystallites !$OMP PARALLEL DO - do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed + do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed myNgrains = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) + + ! --- wind forward --- + if (crystallite_converged(g,i,e)) then if (debugger .and. selectiveDebugger) then !$OMP CRITICAL (write2out) @@ -546,28 +538,33 @@ subroutine crystallite_stressAndItsTangent(updateJaco) endif crystallite_subFrac(g,i,e) = crystallite_subFrac(g,i,e) + crystallite_subStep(g,i,e) formerSubStep = crystallite_subStep(g,i,e) - crystallite_subStep(g,i,e) = min(1.0_pReal-crystallite_subFrac(g,i,e), stepIncreaseCryst*crystallite_subStep(g,i,e)) + crystallite_subStep(g,i,e) = min( 1.0_pReal - crystallite_subFrac(g,i,e), & + stepIncreaseCryst * crystallite_subStep(g,i,e) ) if (crystallite_subStep(g,i,e) > subStepMinCryst) then - crystallite_subTemperature0(g,i,e) = crystallite_Temperature(g,i,e) ! wind forward... - crystallite_subF0(:,:,g,i,e) = crystallite_subF(:,:,g,i,e) ! ...def grad - crystallite_subFp0(:,:,g,i,e) = crystallite_Fp(:,:,g,i,e) ! ...plastic def grad - crystallite_subLp0(:,:,g,i,e) = crystallite_Lp(:,:,g,i,e) ! ...plastic velocity gradient - constitutive_subState0(g,i,e)%p = constitutive_state(g,i,e)%p ! ...microstructure - crystallite_subTstar0_v(:,g,i,e) = crystallite_Tstar_v(:,g,i,e) ! ...2nd PK stress - elseif (formerSubStep > subStepMinCryst) then ! this crystallite just converged + crystallite_subTemperature0(g,i,e) = crystallite_Temperature(g,i,e) ! wind forward... + crystallite_subF0(:,:,g,i,e) = crystallite_subF(:,:,g,i,e) ! ...def grad + crystallite_subFp0(:,:,g,i,e) = crystallite_Fp(:,:,g,i,e) ! ...plastic def grad + crystallite_subLp0(:,:,g,i,e) = crystallite_Lp(:,:,g,i,e) ! ...plastic velocity gradient + constitutive_subState0(g,i,e)%p = constitutive_state(g,i,e)%p ! ...microstructure + crystallite_subTstar0_v(:,g,i,e) = crystallite_Tstar_v(:,g,i,e) ! ...2nd PK stress + elseif (formerSubStep > subStepMinCryst) then ! this crystallite just converged !$OMP CRITICAL (distributionCrystallite) debug_CrystalliteLoopDistribution(min(nCryst+1,NiterationCrystallite)) = & debug_CrystalliteLoopDistribution(min(nCryst+1,NiterationCrystallite)) + 1 !$OMPEND CRITICAL (distributionCrystallite) endif + + ! --- cutback --- + else - crystallite_subStep(g,i,e) = subStepSizeCryst*crystallite_subStep(g,i,e) ! cut step in half and restore... - crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) ! ...temperature - crystallite_Fp(:,:,g,i,e) = crystallite_subFp0(:,:,g,i,e) ! ...plastic def grad - crystallite_invFp(:,:,g,i,e) = math_inv3x3(crystallite_Fp(:,:,g,i,e)) - crystallite_Lp(:,:,g,i,e) = crystallite_subLp0(:,:,g,i,e) ! ...plastic velocity grad - constitutive_state(g,i,e)%p = constitutive_subState0(g,i,e)%p ! ...microstructure - crystallite_Tstar_v(:,g,i,e) = crystallite_subTstar0_v(:,g,i,e) ! ...2nd PK stress + crystallite_subStep(g,i,e) = subStepSizeCryst * crystallite_subStep(g,i,e) ! cut step in half and restore... + crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) ! ...temperature + crystallite_Fp(:,:,g,i,e) = crystallite_subFp0(:,:,g,i,e) ! ...plastic def grad + crystallite_invFp(:,:,g,i,e) = math_inv3x3(crystallite_Fp(:,:,g,i,e)) + crystallite_Lp(:,:,g,i,e) = crystallite_subLp0(:,:,g,i,e) ! ...plastic velocity grad + constitutive_state(g,i,e)%p = constitutive_subState0(g,i,e)%p ! ...microstructure + crystallite_Tstar_v(:,g,i,e) = crystallite_subTstar0_v(:,g,i,e) ! ...2nd PK stress + ! canÕt restore dotState here, since not yet calculated in first cutback after initialization if (debugger .and. selectiveDebugger) then !$OMP CRITICAL (write2out) write(6,'(a78,f10.8)') 'cutback step in crystallite_stressAndItsTangent with new crystallite_subStep: ',& @@ -577,218 +574,53 @@ subroutine crystallite_stressAndItsTangent(updateJaco) endif endif - crystallite_todo(g,i,e) = crystallite_subStep(g,i,e) > subStepMinCryst ! still on track or already done (beyond repair) - if (crystallite_todo(g,i,e)) then ! specify task (according to substep) + ! --- prepare for integration --- + + crystallite_todo(g,i,e) = crystallite_subStep(g,i,e) > subStepMinCryst ! still on track or already done (beyond repair) + if (crystallite_todo(g,i,e)) then crystallite_subF(:,:,g,i,e) = crystallite_subF0(:,:,g,i,e) + & crystallite_subStep(g,i,e) * & (crystallite_partionedF(:,:,g,i,e) - crystallite_partionedF0(:,:,g,i,e)) crystallite_Fe(:,:,g,i,e) = math_mul33x33(crystallite_subF(:,:,g,i,e),crystallite_invFp(:,:,g,i,e)) crystallite_subdt(g,i,e) = crystallite_subStep(g,i,e) * crystallite_dt(g,i,e) - crystallite_converged(g,i,e) = .false. ! start out non-converged + crystallite_converged(g,i,e) = .false. ! start out non-converged endif + enddo enddo enddo !$OMPEND PARALLEL DO - - ! --+>> preguess for state <<+-- - ! - ! incrementing by crystallite_subdt - ! based on constitutive_subState0 - ! results in constitutive_state - ! first loop for collection of state evolution based on old state - ! second loop for updating to new state - - !$OMP PARALLEL DO - do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed - myNgrains = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed - do g = 1,myNgrains - selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) - if (crystallite_todo(g,i,e)) then ! all undone crystallites - call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, & - crystallite_Fp, crystallite_disorientation(:,:,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states - constitutive_previousDotState2(g,i,e)%p = 0.0_pReal - constitutive_previousDotState(g,i,e)%p = 0.0_pReal - constitutive_dotState(g,i,e)%p = 0.0_pReal ! zero out dotStates - endif - enddo; enddo; enddo - !$OMPEND PARALLEL DO - !$OMP PARALLEL DO - do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed - myNgrains = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed - do g = 1,myNgrains - selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) - if (crystallite_todo(g,i,e)) then ! all undone crystallites - call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & - crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & - crystallite_disorientation(:,:,g,i,e), crystallite_subdt(g,i,e), g, i, e) - endif - enddo; enddo; enddo - !$OMPEND PARALLEL DO - crystallite_statedamper = 1.0_pReal - - !$OMP PARALLEL DO - do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed - myNgrains = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed - do g = 1,myNgrains - selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) - if (crystallite_todo(g,i,e)) then ! all undone crystallites - crystallite_stateConverged(g,i,e) = crystallite_updateState(g,i,e) ! update state - crystallite_temperatureConverged(g,i,e) = crystallite_updateTemperature(g,i,e) ! update temperature - if ( .not. crystallite_localConstitution(g,i,e) & - .and. .not. crystallite_todo(g,i,e)) & ! if broken non-local... - crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped - endif - enddo; enddo; enddo - !$OMPEND PARALLEL DO - - ! --+>> state loop <<+-- - - NiterationState = 0_pInt + ! --- integrate --- + + if (any(crystallite_todo)) then + select case(integrator) + case (1) + call crystallite_integrateStateFPI(1) + case (2) + call crystallite_integrateStateEuler(1) + case (3) + call crystallite_integrateStateAdaptiveEuler(1) + case (4) + call crystallite_integrateStateRK4(1) + case(5) + call crystallite_integrateStateRKCK45(1) + endselect + endif - do while ( any(crystallite_todo(:,:,FEsolving_execELem(1):FEsolving_execElem(2))) & - .and. NiterationState < nState) ! convergence loop for crystallite - - NiterationState = NiterationState + 1_pInt - ! --+>> stress integration <<+-- - ! - ! incrementing by crystallite_subdt - ! based on crystallite_subF0,.._subFp0,.._subLp0 - ! constitutive_state is internally interpolated with .._subState0 - ! to account for substepping within _integrateStress - ! results in crystallite_Fp,.._Lp - !$OMP PARALLEL DO - do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed - myNgrains = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed - do g = 1,myNgrains - selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) - if (crystallite_todo(g,i,e)) then ! all undone crystallites - crystallite_todo(g,i,e) = crystallite_integrateStress(1_pInt,g,i,e) ! mode: central solution - if ( .not. crystallite_localConstitution(g,i,e) & - .and. .not. crystallite_todo(g,i,e)) & ! if broken non-local... - crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped - endif - enddo; enddo; enddo - !$OMPEND PARALLEL DO - - if (verboseDebugger) then - !$OMP CRITICAL (write2out) - write(6,*) count(crystallite_todo(:,:,:)),'grains todo after stress integration' - !$OMPEND CRITICAL (write2out) - endif - - ! --+>> state integration <<+-- - ! - ! incrementing by crystallite_subdt - ! based on constitutive_subState0 - ! results in constitutive_state - ! first loop for collection of state evolution based on old state - ! second loop for updating to new state - - !$OMP PARALLEL DO - do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed - myNgrains = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed - do g = 1,myNgrains - if (crystallite_todo(g,i,e)) then ! all undone crystallites - constitutive_previousDotState2(g,i,e)%p = constitutive_previousDotState(g,i,e)%p - constitutive_previousDotState(g,i,e)%p = constitutive_dotState(g,i,e)%p - constitutive_dotState(g,i,e)%p = 0.0_pReal ! zero out dotState - endif - enddo; enddo; enddo - !$OMPEND PARALLEL DO - - crystallite_statedamper = 1.0_pReal - - !$OMP PARALLEL DO - do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed - myNgrains = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed - do g = 1,myNgrains - selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) - if (crystallite_todo(g,i,e)) then ! all undone crystallites - call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & - crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & - crystallite_disorientation(:,:,g,i,e), crystallite_subdt(g,i,e), g, i, e) - dot_prod12 = dot_product( constitutive_dotState(g,i,e)%p - constitutive_previousDotState(g,i,e)%p, & - constitutive_previousDotState(g,i,e)%p - constitutive_previousDotState2(g,i,e)%p ) - dot_prod22 = dot_product( constitutive_previousDotState(g,i,e)%p - constitutive_previousDotState2(g,i,e)%p, & - constitutive_previousDotState(g,i,e)%p - constitutive_previousDotState2(g,i,e)%p ) - if ( dot_prod22 > 0.0_pReal & - .and. ( dot_prod12 < 0.0_pReal & - .or. dot_product(constitutive_dotState(g,i,e)%p, constitutive_previousDotState(g,i,e)%p) < 0.0_pReal) ) & - crystallite_statedamper(g,i,e) = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22) - endif - enddo; enddo; enddo - !$OMPEND PARALLEL DO - - !$OMP PARALLEL DO - do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed - myNgrains = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed - do g = 1,myNgrains - selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) - if (crystallite_todo(g,i,e)) then ! all undone crystallites - crystallite_stateConverged(g,i,e) = crystallite_updateState(g,i,e) ! update state - crystallite_temperatureConverged(g,i,e) = crystallite_updateTemperature(g,i,e) ! update temperature - crystallite_converged(g,i,e) = crystallite_stateConverged(g,i,e) .and. crystallite_temperatureConverged(g,i,e) - if ( .not. crystallite_localConstitution(g,i,e) & - .and. .not. crystallite_todo(g,i,e)) & ! if updateState signals broken non-local... - crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped - if (crystallite_converged(g,i,e)) then - !$OMP CRITICAL (distributionState) - debug_CrystalliteStateLoopDistribution(NiterationState) = & - debug_CrystalliteStateLoopDistribution(NiterationState) + 1 - !$OMPEND CRITICAL (distributionState) - endif - endif - enddo - enddo - enddo - !$OMPEND PARALLEL DO - - if (verboseDebugger) then - !$OMP CRITICAL (write2out) - write(6,*) count(crystallite_converged(:,:,:)),'grains converged after state integration no.', NiterationState - write(6,*) -! write(6,'(8(L,x))') crystallite_converged(:,:,:) -! do e = FEsolving_execElem(1),FEsolving_execElem(2) -! if (any(.not. crystallite_converged(:,:,e))) & -! write(6,'(i4,8(x,L))') e, crystallite_converged(:,:,e) -! enddo - !$OMPEND CRITICAL (write2out) - endif - - if (any(.not. crystallite_converged .and. .not. crystallite_localConstitution)) & ! any non-local not yet converged (or broken)... - crystallite_converged = crystallite_converged .and. crystallite_localConstitution ! ...restart all non-local as not converged - - crystallite_todo = crystallite_todo .and. .not. crystallite_converged ! skip all converged - - if (verboseDebugger) then - !$OMP CRITICAL (write2out) - write(6,*) count(crystallite_converged(:,:,:)),'grains converged after non-local check' - write(6,*) count(crystallite_todo(:,:,:)),'grains todo after state integration no.', NiterationState - write(6,*) - !$OMPEND CRITICAL (write2out) - endif - - enddo ! crystallite convergence loop NiterationCrystallite = NiterationCrystallite + 1 - - enddo ! cutback loop + + enddo ! cutback loop - ! ------ check for non-converged crystallites ------ + + ! --+>> CHECK FOR NON-CONVERGED CRYSTALLITES <<+-- + !$OMP PARALLEL DO - do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed + do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed myNgrains = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains - if (.not. crystallite_converged(g,i,e)) then ! respond fully elastically (might be not required due to becoming terminally ill anyway) + if (.not. crystallite_converged(g,i,e)) then ! respond fully elastically (might be not required due to becoming terminally ill anyway) ! call IO_warning(600,e,i,g) invFp = math_inv3x3(crystallite_partionedFp0(:,:,g,i,e)) Fe_guess = math_mul33x33(crystallite_partionedF(:,:,g,i,e),invFp) @@ -802,11 +634,11 @@ subroutine crystallite_stressAndItsTangent(updateJaco) !$OMPEND PARALLEL DO - ! --+>> stiffness calculation <<+-- + ! --+>> STIFFNESS CALCULATION <<+-- if(updateJaco) then ! Jacobian required - crystallite_statedamper = 1.0_pReal + ! --- BACKUP --- do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed myNgrains = homogenization_Ngrains(mesh_element(3,e)) @@ -827,8 +659,11 @@ subroutine crystallite_stressAndItsTangent(updateJaco) Lp_backup = crystallite_Lp Tstar_v_backup = crystallite_Tstar_v P_backup = crystallite_P - ConvergenceFlag_backup = crystallite_converged + convergenceFlag_backup = crystallite_converged + + ! --- LOCAL STIFFNESS CALCULATION --- + if (all(crystallite_localConstitution) .or. theInc < 1 .or. forceLocalStiffnessCalculation) then ! all grains have local constitution, so local convergence of perturbed grain is sufficient !$OMP PARALLEL DO @@ -836,10 +671,11 @@ subroutine crystallite_stressAndItsTangent(updateJaco) myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains - selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) + selectiveDebugger = .false. ! (e == debug_e .and. i == debug_i .and. g == debug_g) if (crystallite_requested(g,i,e)) then ! first check whether is requested at all! if (crystallite_converged(g,i,e)) then ! grain converged in above iteration - if (debugger .and. selectiveDebugger) then + + if (verboseDebugger .and. selectiveDebugger) then !$OMP CRITICAL (write2out) write (6,*) '#############' write (6,*) 'central solution of cryst_StressAndTangent' @@ -852,57 +688,39 @@ subroutine crystallite_stressAndItsTangent(updateJaco) do perturbation = 1,2 ! forward and backward perturbation if (iand(pert_method,perturbation) > 0) then ! mask for desired direction + dPdF_perturbation(:,:,:,:,perturbation) = crystallite_dPdF0(:,:,:,:,g,i,e) ! initialize stiffness with known good values from last increment myPert = -pert_Fg * (-1.0_pReal)**perturbation ! set perturbation step do k = 1,3; do l = 1,3 ! ...alter individual components crystallite_subF(k,l,g,i,e) = crystallite_subF(k,l,g,i,e) + myPert ! perturb either forward or backward - if (debugger .and. selectiveDebugger) then + + if (verboseDebugger .and. selectiveDebugger) then !$OMP CRITICAL (write2out) write (6,'(a,x,i1,x,i1,x,a)') '[[[[[[[ Stiffness perturbation',k,l,']]]]]]]' write (6,'(a8,3(x,i4),/,3(3(f14.9,x)/))') 'pertF of', g, i, e, crystallite_subF(1:3,:,g,i,e) !$OMPEND CRITICAL (write2out) endif - onTrack = .true. - converged = .false. - NiterationState = 0_pInt - crystallite_statedamper(g,i,e) = 1.0_pReal - do while(.not. converged .and. onTrack .and. NiterationState < nState) ! keep cycling until done (potentially non-converged) - NiterationState = NiterationState + 1_pInt - onTrack = crystallite_integrateStress(2_pInt,g,i,e) ! stress of perturbed situation (overwrites _P,_Tstar_v,_Fp,_Lp,_Fe) with mode:2 - if (onTrack) then - constitutive_dotState(g,i,e)%p = 0.0_pReal - call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & - crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & - crystallite_disorientation(:,:,g,i,e), crystallite_subdt(g,i,e), & - g,i,e) - dot_prod12 = dot_product( & - constitutive_dotState(g,i,e)%p - constitutive_previousDotState(g,i,e)%p, & - constitutive_previousDotState(g,i,e)%p - constitutive_previousDotState2(g,i,e)%p ) - dot_prod22 = dot_product( & - constitutive_previousDotState(g,i,e)%p - constitutive_previousDotState2(g,i,e)%p, & - constitutive_previousDotState(g,i,e)%p - constitutive_previousDotState2(g,i,e)%p ) - if ( dot_prod22 > 0.0_pReal & - .and. ( dot_prod12 < 0.0_pReal & - .or. dot_product(constitutive_dotState(g,i,e)%p, constitutive_previousDotState(g,i,e)%p) & - < 0.0_pReal) ) & - crystallite_statedamper(g,i,e) = 0.75_pReal + & - 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22) - stateConverged = crystallite_updateState(g,i,e) ! update state - temperatureConverged = crystallite_updateTemperature(g,i,e) ! update temperature - converged = stateConverged .and. temperatureConverged - endif - if (debugger .and. selectiveDebugger) then - !$OMP CRITICAL (write2out) - write (6,*) '-------------' - write (6,'(a,x,l,x,l)') 'ontrack + converged:',onTrack,converged - write (6,'(a12,3(x,i4),/,3(3(f12.8,x)/))') 'pertP/MPa of', g, i, e, crystallite_P(1:3,:,g,i,e)/1e6 - write (6,'(a12,3(x,i4),/,3(3(f12.8,x)/))') 'DP/MPa of', g, i, e, & - (crystallite_P(1:3,:,g,i,e)-P_backup(1:3,:,g,i,e))/1e6 - !$OMPEND CRITICAL (write2out) - endif - enddo - if (converged) & ! converged state warrants stiffness update + + ! --- local integration and stiffness calculation --- + + crystallite_converged(g,i,e) = .false. ! start out non-converged + crystallite_todo(g,i,e) = .true. + select case(integratorStiffness) + case (1) + call crystallite_integrateStateFPI(2,g,i,e) + case (2) + call crystallite_integrateStateEuler(2,g,i,e) + case (3) + call crystallite_integrateStateAdaptiveEuler(2,g,i,e) + case (4) + call crystallite_integrateStateRK4(2,g,i,e) + case(5) + call crystallite_integrateStateRKCK45(2,g,i,e) + endselect + if (crystallite_converged(g,i,e)) & ! converged state warrants stiffness update dPdF_perturbation(:,:,k,l,perturbation) = (crystallite_P(:,:,g,i,e) - P_backup(:,:,g,i,e))/myPert ! tangent dP_ij/dFg_kl + ! --- restore --- + mySizeState = constitutive_sizeState(g,i,e) ! number of state variables for this grain mySizeDotState = constitutive_sizeDotState(g,i,e) ! number of dotStates for this grain constitutive_state(g,i,e)%p(1:mySizeState) = constitutive_state_backup(g,i,e)%p(1:mySizeState) @@ -915,10 +733,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) crystallite_Lp(:,:,g,i,e) = Lp_backup(:,:,g,i,e) crystallite_Tstar_v(:,g,i,e) = Tstar_v_backup(:,g,i,e) crystallite_P(:,:,g,i,e) = P_backup(:,:,g,i,e) - !$OMP CRITICAL (out) - debug_StiffnessStateLoopDistribution(NiterationState) = & - debug_StiffnessstateLoopDistribution(NiterationState) + 1 - !$OMPEND CRITICAL (out) + crystallite_converged(g,i,e) = convergenceFlag_backup(g,i,e) + enddo; enddo endif enddo ! perturbation direction @@ -939,91 +755,35 @@ subroutine crystallite_stressAndItsTangent(updateJaco) enddo ! element loop !$OMPEND PARALLEL DO + + ! --- NON-LOCAL STIFFNESS CALCULATION --- + elseif (any(.not. crystallite_localConstitution)) then ! if any nonlocal grain present, we have to do a full loop over all grains after each perturbance + crystallite_dPdF = crystallite_dPdF0 ! initialize stiffness with known good values from last inc + do k = 1,3 do l = 1,3 crystallite_subF(k,l,:,:,:) = crystallite_subF(k,l,:,:,:) + pert_Fg ! perturb single component - NiterationState = 0_pInt + ! --- integration --- + + crystallite_converged = .false. ! start out non-converged crystallite_todo = .true. - do while ( any(crystallite_todo(:,:,FEsolving_execELem(1):FEsolving_execElem(2))) .and. NiterationState < nState) - NiterationState = NiterationState + 1_pInt - - !$OMP PARALLEL DO - do e = FEsolving_execElem(1),FEsolving_execElem(2) - myNgrains = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) - do g = 1,myNgrains - selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) - if (crystallite_todo(g,i,e)) then - crystallite_todo(g,i,e) = crystallite_integrateStress(2_pInt,g,i,e) ! stress integration for stiffness (mode:2) - if ( .not. crystallite_localConstitution(g,i,e) & - .and. .not. crystallite_todo(g,i,e)) & ! if broken non-local... - crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped - endif - enddo; enddo; enddo - !$OMPEND PARALLEL DO - - do e = FEsolving_execElem(1),FEsolving_execElem(2) - myNgrains = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) - do g = 1,myNgrains - if (crystallite_todo(g,i,e)) & - constitutive_dotState(g,i,e)%p = 0.0_pReal ! zero out dotState - enddo; enddo; enddo - - crystallite_statedamper = 1.0_pReal - - !$OMP PARALLEL DO - do e = FEsolving_execElem(1),FEsolving_execElem(2) - myNgrains = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) - do g = 1,myNgrains - selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) - if (crystallite_todo(g,i,e)) then - call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & - crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & - crystallite_disorientation(:,:,g,i,e), crystallite_subdt(g,i,e), & - g,i,e) ! collect dot state - dot_prod12 = dot_product( constitutive_dotState(g,i,e)%p - constitutive_previousDotState(g,i,e)%p, & - constitutive_previousDotState(g,i,e)%p - constitutive_previousDotState2(g,i,e)%p ) - dot_prod22 = dot_product( constitutive_previousDotState(g,i,e)%p - constitutive_previousDotState2(g,i,e)%p,& - constitutive_previousDotState(g,i,e)%p - constitutive_previousDotState2(g,i,e)%p ) - if ( dot_prod22 > 0.0_pReal & - .and. ( dot_prod12 < 0.0_pReal & - .or. dot_product(constitutive_dotState(g,i,e)%p, & - constitutive_previousDotState(g,i,e)%p) < 0.0_pReal) ) & - crystallite_statedamper(g,i,e) = 0.75_pReal & - + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22) - endif - enddo; enddo; enddo - !$OMPEND PARALLEL DO - - !$OMP PARALLEL DO - do e = FEsolving_execElem(1),FEsolving_execElem(2) - myNgrains = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) - do g = 1,myNgrains - selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) - if (crystallite_todo(g,i,e)) then - crystallite_stateConverged(g,i,e) = crystallite_updateState(g,i,e) ! update state - crystallite_temperatureConverged(g,i,e) = crystallite_updateTemperature(g,i,e) ! update temperature - crystallite_converged(g,i,e) = crystallite_stateConverged(g,i,e) & - .and. crystallite_temperatureConverged(g,i,e) - if ( .not. crystallite_localConstitution(g,i,e) & - .and. .not. crystallite_todo(g,i,e)) & ! if updateState signals broken non-local... - crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped - endif - enddo; enddo; enddo - !$OMPEND PARALLEL DO - - if (any(.not. crystallite_converged .and. .not. crystallite_localConstitution)) & ! any non-local not yet converged? - crystallite_converged = crystallite_converged .and. crystallite_localConstitution ! all non-local not converged - - crystallite_todo = crystallite_todo .and. .not. crystallite_converged ! skip all converged - - enddo + select case(integratorStiffness) + case (1) + call crystallite_integrateStateFPI(2) + case (2) + call crystallite_integrateStateEuler(2) + case (3) + call crystallite_integrateStateAdaptiveEuler(2) + case (4) + call crystallite_integrateStateRK4(2) + case(5) + call crystallite_integrateStateRKCK45(2) + endselect + + ! --- stiffness calculation --- do e = FEsolving_execElem(1),FEsolving_execElem(2) myNgrains = homogenization_Ngrains(mesh_element(3,e)) @@ -1031,15 +791,13 @@ subroutine crystallite_stressAndItsTangent(updateJaco) do g = 1,myNgrains if (crystallite_converged(g,i,e)) then ! if stiffness calculation converged... crystallite_dPdF(:,:,k,l,g,i,e) = (crystallite_P(:,:,g,i,e) - P_backup(:,:,g,i,e))/pert_Fg ! ... use tangent dP_ij/dFg_kl - !$OMP CRITICAL (out) - debug_StiffnessStateLoopDistribution(NiterationState) = & - debug_StiffnessstateLoopDistribution(NiterationState) + 1 - !$OMPEND CRITICAL (out) - elseif (.not. ConvergenceFlag_backup(g,i,e)) then ! if crystallite didnÕt converge before... + elseif (.not. convergenceFlag_backup(g,i,e)) then ! if crystallite didnt converge before... crystallite_dPdF(:,:,:,:,g,i,e) = crystallite_fallbackdPdF(:,:,:,:,g,i,e) ! ... use (elastic) fallback endif enddo; enddo; enddo - + + ! --- restore --- + do e = FEsolving_execElem(1),FEsolving_execElem(2) myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) @@ -1057,8 +815,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) crystallite_Lp = Lp_backup crystallite_Tstar_v = Tstar_v_backup crystallite_P = P_backup - - crystallite_converged = ConvergenceFlag_backup + crystallite_converged = convergenceFlag_backup enddo;enddo ! k,l loop @@ -1070,98 +827,1416 @@ endsubroutine +!******************************************************************** +! integrate stress, state and Temperature with +! 4h order explicit Runge Kutta method +!******************************************************************** +subroutine crystallite_integrateStateRK4(mode,gg,ii,ee) + +!*** variables and functions from other modules ***! +use prec, only: pInt, & + pReal +use debug, only: debugger, & + selectiveDebugger, & + verboseDebugger, & + debug_e, & + debug_i, & + debug_g, & + debug_StateLoopDistribution +use FEsolving, only: FEsolving_execElem, & + FEsolving_execIP +use mesh, only: mesh_element, & + mesh_NcpElems, & + mesh_maxNips +use material, only: homogenization_Ngrains, & + homogenization_maxNgrains +use constitutive, only: constitutive_sizeDotState, & + constitutive_state, & + constitutive_subState0, & + constitutive_dotState, & + constitutive_RK4dotState, & + constitutive_collectDotState, & + constitutive_dotTemperature, & + constitutive_microstructure + +implicit none + +real(pReal), dimension(4), parameter :: timeStepFraction = (/0.5_pReal, 0.5_pReal, 1.0_pReal, 1.0_pReal/) ! weight of slope used for Runge Kutta integration +real(pReal), dimension(3), parameter :: weight = (/2.0_pReal, 2.0_pReal, 1.0_pReal/) ! factor giving the fraction of the original timestep used for Runge Kutta Integration + +!*** input variables ***! +integer(pInt), intent(in) :: mode ! mode of calculation; 1: central solution, 2: stiffness (by perturbation) +integer(pInt), optional, intent(in):: ee, & ! element index + ii, & ! integration point index + gg ! grain index + +!*** output variables ***! + +!*** local variables ***! +integer(pInt) e, & ! element index in element loop + i, & ! integration point index in ip loop + g, & ! grain index in grain loop + n, & + mySizeDotState +integer(pInt), dimension(2) :: eIter ! bounds for element iteration +integer(pInt), dimension(2,mesh_NcpElems) :: iIter, & ! bounds for ip iteration + gIter ! bounds for grain iteration +real(pReal), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & + RK4dotTemperature ! evolution of Temperature of each grain for Runge Kutta integration +logical singleRun ! flag indicating computation for single (g,i,e) triple + + +if (present(ee) .and. present(ii) .and. present(gg)) then + eIter = ee + iIter(:,ee) = ii + gIter(:,ee) = gg + singleRun = .true. +else + eIter = FEsolving_execElem(1:2) + do e = eIter(1),eIter(2) + iIter(:,e) = FEsolving_execIP(1:2,e) + gIter(:,e) = (/1,homogenization_Ngrains(mesh_element(3,e))/) + enddo + singleRun = .false. +endif + + +! --- UPDATE DEPENDENT STATES --- + +!$OMP PARALLEL DO + do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, & + crystallite_Fp, crystallite_disorientation(:,:,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states + endif + enddo; enddo; enddo +!$OMPEND PARALLEL DO + + +! --- FIRST RUNGE KUTTA STEP --- + +RK4dotTemperature = 0.0_pReal +!$OMP PARALLEL DO + do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + selectiveDebugger = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g) + call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & + crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & + crystallite_disorientation(:,:,g,i,e), crystallite_subdt(g,i,e), g,i,e) + crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), & + crystallite_Temperature(g,i,e),g,i,e) + if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState + .or. crystallite_dotTemperature(g,i,e)/=crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature + if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local... + !$OMP CRITICAL + crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped + !$OMPEND CRITICAL + else ! if broken local... + crystallite_todo(g,i,e) = .false. ! ... skip this one next time + endif + else ! everything is fine + constitutive_RK4dotState(g,i,e)%p = constitutive_dotState(g,i,e)%p ! initial contribution to RK slope + RK4dotTemperature(g,i,e) = crystallite_dotTemperature(g,i,e) + endif + endif + enddo; enddo; enddo +!$OMPEND PARALLEL DO + + +! --- SECOND TO FOURTH RUNGE KUTTA STEP PLUS FINAL INTEGRATION --- + +do n = 1,4 + + + ! --- state update --- + + !$OMP PARALLEL DO + do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + mySizeDotState = constitutive_sizeDotState(g,i,e) + constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_subState0(g,i,e)%p(1:mySizeDotState) & + + constitutive_dotState(g,i,e)%p(1:mySizeDotState) * crystallite_subdt(g,i,e) * timeStepFraction(n) + crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) & + + crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e) * timeStepFraction(n) + endif + enddo; enddo; enddo + !$OMPEND PARALLEL DO + + + ! --- update dependent states --- + + !$OMP PARALLEL DO + do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, & + crystallite_Fp, crystallite_disorientation(:,:,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states + endif + enddo; enddo; enddo + !$OMPEND PARALLEL DO + + + ! --- stress integration --- + + !$OMP PARALLEL DO + do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + selectiveDebugger = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g) + if (crystallite_integrateStress(mode,g,i,e,timeStepFraction(n))) then ! fraction of original times step + if (n == 4) then ! final integration step + if (mode==1 .and. verboseDebugger .and. e == debug_e .and. i == debug_i .and. g == debug_g) then + mySizeDotState = constitutive_sizeDotState(g,i,e) + !$OMP CRITICAL (write2out) + write(6,*) '::: updateState',g,i,e + write(6,*) + write(6,'(a,/,12(e14.8,x))') 'updateState: dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState) + write(6,*) + write(6,'(a,/,12(e14.8,x))') 'updateState: new state', constitutive_state(g,i,e)%p(1:mySizeDotState) + write(6,*) + !$OMPEND CRITICAL (write2out) + endif + crystallite_converged(g,i,e) = .true. ! ... converged per definition + crystallite_todo(g,i,e) = .false. ! ... integration done + !$OMP CRITICAL (distributionState) + debug_StateLoopDistribution(n,mode) = debug_StateLoopDistribution(n,mode) + 1 + !$OMPEND CRITICAL (distributionState) + endif + else ! broken stress integration + if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local... + !$OMP CRITICAL + crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped + !$OMPEND CRITICAL + else ! if broken local... + crystallite_todo(g,i,e) = .false. ! ... skip this one next time + endif + endif + endif + enddo; enddo; enddo + !$OMPEND PARALLEL DO + + + ! --- dot state and RK dot state--- + + if (n < 4) then + !$OMP PARALLEL DO + do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + selectiveDebugger = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g) + call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & + crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & + crystallite_disorientation(:,:,g,i,e), & + timeStepFraction(n)*crystallite_subdt(g,i,e), g,i,e) ! fraction of original timestep + crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), & + crystallite_Temperature(g,i,e),g,i,e) + if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState + .or. crystallite_dotTemperature(g,i,e)/=crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature + if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local... + !$OMP CRITICAL + crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped + !$OMPEND CRITICAL + else ! if broken local... + crystallite_todo(g,i,e) = .false. ! ... skip this one next time + endif + else ! everything is fine + constitutive_RK4dotState(g,i,e)%p = constitutive_RK4dotState(g,i,e)%p + weight(n)*constitutive_dotState(g,i,e)%p + RK4dotTemperature(g,i,e) = RK4dotTemperature(g,i,e) + weight(n)*crystallite_dotTemperature(g,i,e) + if (n == 3) then + constitutive_dotState(g,i,e)%p = constitutive_RK4dotState(g,i,e)%p / 6.0_pReal ! use weighted RKdotState for final integration + endif + endif + endif + enddo; enddo; enddo + !$OMPEND PARALLEL DO + endif + +enddo + + +! --- CHECK CONVERGENCE --- + +crystallite_todo = .false. ! done with integration +if ( .not. (mode == 2 .and. singleRun) & ! except for local stiffness calculation: + .and. any(.not. crystallite_converged .and. .not. crystallite_localConstitution)) then ! any non-local not yet converged (or broken)... + crystallite_converged = crystallite_converged .and. crystallite_localConstitution ! ...restart all non-local as not converged +endif + +endsubroutine + + + +!******************************************************************** +! integrate stress, state and Temperature with +! 5th order Runge-Kutta Cash-Karp method with adaptive step size +! (use 5th order solution to advance = "local extrapolation") +!******************************************************************** +subroutine crystallite_integrateStateRKCK45(mode,gg,ii,ee) + +!*** variables and functions from other modules ***! +use prec, only: pInt, & + pReal +use debug, only: debugger, & + selectiveDebugger, & + verboseDebugger, & + debug_e, & + debug_i, & + debug_g, & + debug_StateLoopDistribution +use numerics, only: rTol_crystalliteState, & + rTol_crystalliteTemperature, & + subStepSizeCryst, & + stepIncreaseCryst +use FEsolving, only: FEsolving_execElem, & + FEsolving_execIP +use mesh, only: mesh_element, & + mesh_NcpElems, & + mesh_maxNips +use material, only: homogenization_Ngrains, & + homogenization_maxNgrains +use constitutive, only: constitutive_sizeDotState, & + constitutive_maxSizeDotState, & + constitutive_state, & + constitutive_relevantState, & + constitutive_subState0, & + constitutive_dotState, & + constitutive_RKCK45dotState, & + constitutive_collectDotState, & + constitutive_dotTemperature, & + constitutive_microstructure + +implicit none + + +!*** input variables ***! +integer(pInt), intent(in) :: mode ! mode of calculation; 1: central solution, 2: stiffness (by perturbation) +integer(pInt), optional, intent(in):: ee, & ! element index + ii, & ! integration point index + gg ! grain index + +!*** output variables ***! + +!*** local variables ***! +integer(pInt) e, & ! element index in element loop + i, & ! integration point index in ip loop + g, & ! grain index in grain loop + j, & + n, & ! stage index in integration stage loop + sizeDotState, & ! size of dot State + s ! state index +integer(pInt), dimension(2) :: eIter ! bounds for element iteration +integer(pInt), dimension(2,mesh_NcpElems) :: iIter, & ! bounds for ip iteration + gIter ! bounds for grain iteration +real(pReal), dimension(6,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & + RKCK45dotTemperature ! evolution of Temperature of each grain for Runge Kutta Cash Karp integration +real(pReal), dimension(5,5) :: a ! coefficients in Butcher tableau (used for preliminary integration in stages 2 to 6) +real(pReal), dimension(6) :: b, db ! coefficients in Butcher tableau (used for final integration and error estimate) +real(pReal), dimension(5) :: c ! coefficients in Butcher tableau (fractions of original time step in stages 2 to 6) +real(pReal), dimension(constitutive_maxSizeDotState,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & + stateResiduum, & ! residuum from evolution in micrstructure + relStateResiduum ! relative residuum from evolution in microstructure +real(pReal), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & + temperatureResiduum, & ! residuum from evolution in temperature + relTemperatureResiduum ! relative residuum from evolution in temperature +logical singleRun ! flag indicating computation for single (g,i,e) triple + + +! --- FILL BUTCHER TABLEAU --- + +a = 0.0_pReal +b = 0.0_pReal +db = 0.0_pReal +c = 0.0_pReal + +a(1,1) = 0.2_pReal +a(1,2) = 0.075_pReal +a(2,2) = 0.225_pReal +a(1,3) = 0.3_pReal +a(2,3) = -0.9_pReal +a(3,3) = 1.2_pReal +a(1,4) = -11.0_pReal / 54.0_pReal +a(2,4) = 2.5_pReal +a(3,4) = -70.0_pReal / 27.0_pReal +a(4,4) = 35.0_pReal / 27.0_pReal +a(1,5) = 1631.0_pReal / 55296.0_pReal +a(2,5) = 175.0_pReal / 512.0_pReal +a(3,5) = 575.0_pReal / 13824.0_pReal +a(4,5) = 44275.0_pReal / 110592.0_pReal +a(5,5) = 253.0_pReal / 4096.0_pReal + +b(1) = 37.0_pReal / 378.0_pReal +b(3) = 250.0_pReal / 621.0_pReal +b(4) = 125.0_pReal / 594.0_pReal +b(6) = 512.0_pReal / 1771.0_pReal + +db(1) = b(1) - 2825.0_pReal / 27648.0_pReal +db(3) = b(3) - 18575.0_pReal / 48384.0_pReal +db(4) = b(4) - 13525.0_pReal / 55296.0_pReal +db(5) = - 277.0_pReal / 14336.0_pReal +db(6) = b(6) - 0.25_pReal + +c(1) = 0.2_pReal +c(2) = 0.3_pReal +c(3) = 0.6_pReal +c(4) = 1.0_pReal +c(5) = 0.875_pReal + + +! --- LOOP ITERATOR FOR ELEMENT, GRAIN, IP --- + +if (present(ee) .and. present(ii) .and. present(gg)) then + eIter = ee + iIter(:,ee) = ii + gIter(:,ee) = gg + singleRun = .true. +else + eIter = FEsolving_execElem(1:2) + do e = eIter(1),eIter(2) + iIter(:,e) = FEsolving_execIP(1:2,e) + gIter(:,e) = (/1,homogenization_Ngrains(mesh_element(3,e))/) + enddo + singleRun = .false. +endif + + +! --- UPDATE DEPENDENT STATES --- + +!$OMP PARALLEL DO + do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, & + crystallite_Fp, crystallite_disorientation(:,:,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states + endif + enddo; enddo; enddo +!$OMPEND PARALLEL DO + + +! --- FIRST RUNGE KUTTA STEP --- + +!$OMP PARALLEL DO + do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + selectiveDebugger = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g) + call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & + crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & + crystallite_disorientation(:,:,g,i,e), crystallite_subdt(g,i,e), g,i,e) + crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), & + crystallite_Temperature(g,i,e),g,i,e) + if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState + .or. crystallite_dotTemperature(g,i,e)/=crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature + if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local... + !$OMP CRITICAL + crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped + !$OMPEND CRITICAL + else ! if broken local... + crystallite_todo(g,i,e) = .false. ! ... skip this one next time + endif + else ! everything is fine + constitutive_RKCK45dotState(1,g,i,e)%p = constitutive_dotState(g,i,e)%p ! initial contribution to RK slope + RKCK45dotTemperature(1,g,i,e) = crystallite_dotTemperature(g,i,e) + endif + endif + enddo; enddo; enddo +!$OMPEND PARALLEL DO + + +! --- SECOND TO SIXTH RUNGE KUTTA STEP --- + +do n = 1,5 + + + ! --- state update --- + + !$OMP PARALLEL DO + do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + sizeDotState = constitutive_sizeDotState(g,i,e) + constitutive_dotState(g,i,e)%p = 0.0_pReal + crystallite_dotTemperature(g,i,e) = 0.0_pReal + do j = 1,n + constitutive_dotState(g,i,e)%p = constitutive_dotState(g,i,e)%p + a(j,n) * constitutive_RKCK45dotState(j,g,i,e)%p + crystallite_dotTemperature(g,i,e) = crystallite_dotTemperature(g,i,e) + a(j,n) * RKCK45dotTemperature(j,g,i,e) + enddo + constitutive_state(g,i,e)%p(1:sizeDotState) = constitutive_subState0(g,i,e)%p(1:sizeDotState) & + + constitutive_dotState(g,i,e)%p(1:sizeDotState) * crystallite_subdt(g,i,e) + crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) & + + crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e) + endif + enddo; enddo; enddo + !$OMPEND PARALLEL DO + + + ! --- update dependent states --- + + !$OMP PARALLEL DO + do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, & + crystallite_Fp, crystallite_disorientation(:,:,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states + endif + enddo; enddo; enddo + !$OMPEND PARALLEL DO + + + ! --- stress integration --- + + !$OMP PARALLEL DO + do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + selectiveDebugger = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g .and. mode == 1) + if (crystallite_todo(g,i,e)) then + if (.not. crystallite_integrateStress(mode,g,i,e,c(n))) then ! fraction of original time step + if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local... + !$OMP CRITICAL + crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped + !$OMPEND CRITICAL + else ! if broken local... + crystallite_todo(g,i,e) = .false. ! ... skip this one next time + endif + endif + endif + enddo; enddo; enddo + !$OMPEND PARALLEL DO + + + ! --- dot state and RK dot state--- + + !$OMP PARALLEL DO + do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + selectiveDebugger = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g) + if (crystallite_todo(g,i,e)) then + call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & + crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & + crystallite_disorientation(:,:,g,i,e), c(n)*crystallite_subdt(g,i,e), g,i,e) ! fraction of original timestep + crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), & + crystallite_Temperature(g,i,e),g,i,e) + if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState + .or. crystallite_dotTemperature(g,i,e)/=crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature + if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local... + !$OMP CRITICAL + crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped + !$OMPEND CRITICAL + else ! if broken local... + crystallite_todo(g,i,e) = .false. ! ... skip this one next time + endif + else ! everything is fine + constitutive_RKCK45dotState(n+1,g,i,e)%p = constitutive_dotState(g,i,e)%p + RKCK45dotTemperature(n+1,g,i,e) = crystallite_dotTemperature(g,i,e) + endif + endif + enddo; enddo; enddo + !$OMPEND PARALLEL DO + +enddo + + +! --- STATE UPDATE WITH ERROR ESTIMATE FOR STATE AND TEMPERATURE --- + +stateResiduum = 0.0_pReal +temperatureResiduum = 0.0_pReal +relStateResiduum = 0.0_pReal +relTemperatureResiduum = 0.0_pReal +!$OMP PARALLEL DO + do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode==1) + sizeDotState = constitutive_sizeDotState(g,i,e) + + + ! --- absolute residuum in state and temperature --- + + do j = 1,6 + stateResiduum(1:sizeDotState,g,i,e) = stateResiduum(1:sizeDotState,g,i,e) & + + db(j) * constitutive_RKCK45dotState(j,g,i,e)%p(1:sizeDotState) * crystallite_subdt(g,i,e) + temperatureResiduum(g,i,e) = temperatureResiduum(g,i,e) + db(j) * RKCK45dotTemperature(j,g,i,e) * crystallite_subdt(g,i,e) + enddo + + + ! --- dot state and dot temperature --- + + constitutive_dotState(g,i,e)%p = 0.0_pReal + crystallite_dotTemperature(g,i,e) = 0.0_pReal + do j = 1,6 + constitutive_dotState(g,i,e)%p = constitutive_dotState(g,i,e)%p + b(j) * constitutive_RKCK45dotState(j,g,i,e)%p + crystallite_dotTemperature(g,i,e) = crystallite_dotTemperature(g,i,e) + b(j) * RKCK45dotTemperature(j,g,i,e) + enddo + + + ! --- state and temperature update and relative residui --- + + constitutive_state(g,i,e)%p(1:sizeDotState) = constitutive_subState0(g,i,e)%p(1:sizeDotState) & + + constitutive_dotState(g,i,e)%p(1:sizeDotState) * crystallite_subdt(g,i,e) + crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) & + + crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e) + forall (s = 1:sizeDotState, constitutive_state(g,i,e)%p(s) > constitutive_relevantState(g,i,e)%p(s)) & + relStateResiduum(s,g,i,e) = abs(stateResiduum(s,g,i,e)) / constitutive_state(g,i,e)%p(s) / rTol_crystalliteState + if (crystallite_Temperature(g,i,e) > 0) & + relTemperatureResiduum(g,i,e) = abs(temperatureResiduum(g,i,e)) / crystallite_Temperature(g,i,e) & + / rTol_crystalliteTemperature + + if (verboseDebugger .and. selectiveDebugger) then + !$OMP CRITICAL (write2out) + write(6,*) '::: updateState',g,i,e + write(6,*) + write(6,'(a,/,12(f12.1,x))') 'updateState: absolute residuum', stateResiduum(1:sizeDotState,g,i,e) + write(6,*) + write(6,'(a,/,12(f12.1,x))') 'updateState: resid tolerance', relStateResiduum(1:sizeDotState,g,i,e) + write(6,*) +! write(6,'(a)') 'updateState: RKCK45dotState' +! do j = 1,6 +! write(6,'(12(e14.8,x))') constitutive_RKCK45dotState(j,g,i,e)%p(1:sizeDotState) +! write(6,*) +! enddo + write(6,'(a,/,12(e14.8,x))') 'updateState: dotState', constitutive_dotState(g,i,e)%p(1:sizeDotState) + write(6,*) + write(6,'(a,/,12(e14.8,x))') 'updateState: new state', constitutive_state(g,i,e)%p(1:sizeDotState) + write(6,*) + !$OMPEND CRITICAL (write2out) + endif + + endif + enddo; enddo; enddo +!$OMPEND PARALLEL DO + + +! --- UPDATE DEPENDENT STATES IF RESIDUUM BELOW TOLERANCE --- + +!$OMP PARALLEL DO + do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + sizeDotState = constitutive_sizeDotState(g,i,e) + if ( all(relStateResiduum(1:sizeDotState,g,i,e) < 1.0_pReal) .and. relTemperatureResiduum(g,i,e) < 1.0_pReal ) then + call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, & + crystallite_Fp, crystallite_disorientation(:,:,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states + endif + endif + enddo; enddo; enddo +!$OMPEND PARALLEL DO + + +! --- FINAL STRESS INTEGRATION STEP IF RESIDUUM BELOW TOLERANCE --- + +!$OMP PARALLEL DO + do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + selectiveDebugger = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g .and. mode==1) + sizeDotState = constitutive_sizeDotState(g,i,e) + if ( all(relStateResiduum(1:sizeDotState,g,i,e) < 1.0_pReal) .and. relTemperatureResiduum(g,i,e) < 1.0_pReal ) then + + if (crystallite_integrateStress(mode,g,i,e)) then + crystallite_converged(g,i,e) = .true. ! ... converged per definitionem + crystallite_todo(g,i,e) = .false. ! ... integration done + !$OMP CRITICAL (distributionState) + debug_StateLoopDistribution(6,mode) = debug_StateLoopDistribution(6,mode) + 1 + !$OMPEND CRITICAL (distributionState) + else + if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local... + !$OMP CRITICAL + crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped + !$OMPEND CRITICAL + endif + endif + + endif + endif + enddo; enddo; enddo +!$OMPEND PARALLEL DO + + +! --- nonlocal convergence check --- + +if (verboseDebugger .and. mode==1) write(6,*) 'crystallite_converged',crystallite_converged +if ( .not. (mode == 2 .and. singleRun) ) then ! except for local stiffness calculation: + + if ( any(.not. crystallite_converged .and. .not. crystallite_localConstitution)) then ! any non-local not yet converged (or broken)... + crystallite_converged = crystallite_converged .and. crystallite_localConstitution ! ...restart all non-local as not converged + endif + +endif + +endsubroutine + + + +!******************************************************************** +! integrate stress, state and Temperature with +! 1nd order Euler method with adaptive step size +!******************************************************************** +subroutine crystallite_integrateStateAdaptiveEuler(mode,gg,ii,ee) + +!*** variables and functions from other modules ***! +use prec, only: pInt, & + pReal +use debug, only: debugger, & + selectiveDebugger, & + verboseDebugger, & + debug_e, & + debug_i, & + debug_g, & + debug_StateLoopDistribution +use numerics, only: rTol_crystalliteState, & + rTol_crystalliteTemperature, & + subStepSizeCryst, & + stepIncreaseCryst +use FEsolving, only: FEsolving_execElem, & + FEsolving_execIP +use mesh, only: mesh_element, & + mesh_NcpElems, & + mesh_maxNips +use material, only: homogenization_Ngrains, & + homogenization_maxNgrains +use constitutive, only: constitutive_sizeDotState, & + constitutive_maxSizeDotState, & + constitutive_state, & + constitutive_relevantState, & + constitutive_subState0, & + constitutive_dotState, & + constitutive_collectDotState, & + constitutive_dotTemperature, & + constitutive_microstructure + +implicit none + + +!*** input variables ***! +integer(pInt), intent(in) :: mode ! mode of calculation; 1: central solution, 2: stiffness (by perturbation) +integer(pInt), optional, intent(in):: ee, & ! element index + ii, & ! integration point index + gg ! grain index + +!*** output variables ***! + +!*** local variables ***! +integer(pInt) e, & ! element index in element loop + i, & ! integration point index in ip loop + g, & ! grain index in grain loop + j, & + n, & ! stage index in integration stage loop + sizeDotState, & ! size of dot State + s ! state index +integer(pInt), dimension(2) :: eIter ! bounds for element iteration +integer(pInt), dimension(2,mesh_NcpElems) :: iIter, & ! bounds for ip iteration + gIter ! bounds for grain iteration +real(pReal), dimension(constitutive_maxSizeDotState,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & + stateResiduum, & ! residuum from evolution in micrstructure + relStateResiduum ! relative residuum from evolution in microstructure +real(pReal), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & + temperatureResiduum, & ! residuum from evolution in temperature + relTemperatureResiduum ! relative residuum from evolution in temperature +logical singleRun ! flag indicating computation for single (g,i,e) triple + + +! --- LOOP ITERATOR FOR ELEMENT, GRAIN, IP --- + +if (present(ee) .and. present(ii) .and. present(gg)) then + eIter = ee + iIter(:,ee) = ii + gIter(:,ee) = gg + singleRun = .true. +else + eIter = FEsolving_execElem(1:2) + do e = eIter(1),eIter(2) + iIter(:,e) = FEsolving_execIP(1:2,e) + gIter(:,e) = (/1,homogenization_Ngrains(mesh_element(3,e))/) + enddo + singleRun = .false. +endif + + +! --- UPDATE DEPENDENT STATES (EULER INTEGRATION) --- + +!$OMP PARALLEL DO + do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, & + crystallite_Fp, crystallite_disorientation(:,:,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states + endif + enddo; enddo; enddo +!$OMPEND PARALLEL DO + + +! --- DOT STATE AND TEMPERATURE (EULER INTEGRATION) --- + +!$OMP PARALLEL DO + do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + selectiveDebugger = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g) + + call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & + crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & + crystallite_disorientation(:,:,g,i,e), crystallite_subdt(g,i,e), g,i,e) + crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), & + crystallite_Temperature(g,i,e),g,i,e) + + if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState + .or. crystallite_dotTemperature(g,i,e)/=crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature + if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local... + !$OMP CRITICAL + crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped + !$OMPEND CRITICAL + else ! if broken local... + crystallite_todo(g,i,e) = .false. ! ... skip this one next time + endif + else + stateResiduum(:,g,i,e) = - 0.5_pReal * constitutive_dotState(g,i,e)%p * crystallite_subdt(g,i,e) ! contribution to absolute residuum in state and temperature + temperatureResiduum(g,i,e) = - 0.5_pReal * crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e) + endif + endif + enddo; enddo; enddo +!$OMPEND PARALLEL DO + + +! --- STATE UPDATE (EULER INTEGRATION) --- + +!$OMP PARALLEL DO + do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + sizeDotState = constitutive_sizeDotState(g,i,e) + constitutive_state(g,i,e)%p(1:sizeDotState) = constitutive_subState0(g,i,e)%p(1:sizeDotState) & + + constitutive_dotState(g,i,e)%p(1:sizeDotState) * crystallite_subdt(g,i,e) + crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) & + + crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e) + endif + enddo; enddo; enddo +!$OMPEND PARALLEL DO + + +! --- UPDATE DEPENDENT STATES (EULER INTEGRATION) --- + +!$OMP PARALLEL DO + do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, & + crystallite_Fp, crystallite_disorientation(:,:,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states + endif + enddo; enddo; enddo +!$OMPEND PARALLEL DO + + +! --- STRESS INTEGRATION (EULER INTEGRATION) --- + +!$OMP PARALLEL DO + do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + selectiveDebugger = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g .and. mode == 1) + if (crystallite_todo(g,i,e)) then + if (.not. crystallite_integrateStress(mode,g,i,e)) then + if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local... + !$OMP CRITICAL + crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped + !$OMPEND CRITICAL + else ! if broken local... + crystallite_todo(g,i,e) = .false. ! ... skip this one next time + endif + endif + endif + enddo; enddo; enddo +!$OMPEND PARALLEL DO + + +! --- DOT STATE AND TEMPERATURE (HEUN METHOD) --- + +!$OMP PARALLEL DO + do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + selectiveDebugger = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g) + if (crystallite_todo(g,i,e)) then + call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & + crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & + crystallite_disorientation(:,:,g,i,e), crystallite_subdt(g,i,e), g,i,e) + crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), & + crystallite_Temperature(g,i,e),g,i,e) + if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState + .or. crystallite_dotTemperature(g,i,e)/=crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature + if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local... + !$OMP CRITICAL + crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped + !$OMPEND CRITICAL + else ! if broken local... + crystallite_todo(g,i,e) = .false. ! ... skip this one next time + endif + endif + endif + enddo; enddo; enddo +!$OMPEND PARALLEL DO + + +! --- ERROR ESTIMATE FOR STATE AND TEMPERATURE (HEUN METHOD) --- + +relStateResiduum = 0.0_pReal +relTemperatureResiduum = 0.0_pReal +!$OMP PARALLEL DO + do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode==1) + sizeDotState = constitutive_sizeDotState(g,i,e) + + + ! --- contribution of heun step to absolute residui --- + + stateResiduum(:,g,i,e) = stateResiduum(:,g,i,e) & + + 0.5_pReal * constitutive_dotState(g,i,e)%p * crystallite_subdt(g,i,e) ! contribution to absolute residuum in state and temperature + temperatureResiduum(g,i,e) = temperatureResiduum(g,i,e) & + + 0.5_pReal * crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e) + + + ! --- relative residui --- + + forall (s = 1:sizeDotState, constitutive_state(g,i,e)%p(s) > constitutive_relevantState(g,i,e)%p(s)) & + relStateResiduum(s,g,i,e) = abs(stateResiduum(s,g,i,e)) / constitutive_state(g,i,e)%p(s) / rTol_crystalliteState + if (crystallite_Temperature(g,i,e) > 0) & + relTemperatureResiduum(g,i,e) = abs(temperatureResiduum(g,i,e)) / crystallite_Temperature(g,i,e) & + / rTol_crystalliteTemperature + + if (verboseDebugger .and. selectiveDebugger) then + !$OMP CRITICAL (write2out) + write(6,*) '::: updateState',g,i,e + write(6,*) + write(6,'(a,/,12(f12.1,x))') 'updateState: absolute residuum', stateResiduum(1:sizeDotState,g,i,e) + write(6,*) + write(6,'(a,/,12(f12.1,x))') 'updateState: resid tolerance', relStateResiduum(1:sizeDotState,g,i,e) + write(6,*) + write(6,'(a,/,12(e14.8,x))') 'updateState: dotState', constitutive_dotState(g,i,e)%p(1:sizeDotState) & + - 2.0_pReal * stateResiduum(1:sizeDotState,g,i,e) / crystallite_subdt(g,i,e) ! calculate former dotstate from higher order solution and state residuum + write(6,*) + write(6,'(a,/,12(e14.8,x))') 'updateState: new state', constitutive_state(g,i,e)%p(1:sizeDotState) + write(6,*) + !$OMPEND CRITICAL (write2out) + endif + + + ! --- converged ? --- + + if ( all(relStateResiduum(1:sizeDotState,g,i,e) < 1.0_pReal) .and. relTemperatureResiduum(g,i,e) < 1.0_pReal ) then + + crystallite_converged(g,i,e) = .true. ! ... converged per definitionem + crystallite_todo(g,i,e) = .false. ! ... integration done + !$OMP CRITICAL (distributionState) + debug_StateLoopDistribution(6,mode) = debug_StateLoopDistribution(6,mode) + 1 + !$OMPEND CRITICAL (distributionState) + + endif + + endif + enddo; enddo; enddo +!$OMPEND PARALLEL DO + + +! --- NONLOCAL CONVERGENCE CHECK --- + +if (verboseDebugger .and. mode==1) write(6,*) 'crystallite_converged',crystallite_converged +if ( .not. (mode == 2 .and. singleRun) ) then ! except for local stiffness calculation: + if ( any(.not. crystallite_converged .and. .not. crystallite_localConstitution)) then ! any non-local not yet converged (or broken)... + crystallite_converged = crystallite_converged .and. crystallite_localConstitution ! ...restart all non-local as not converged + endif +endif + +endsubroutine + + + +!******************************************************************** +! integrate stress, state and Temperature with +! 1st order explicit Euler method +!******************************************************************** +subroutine crystallite_integrateStateEuler(mode,gg,ii,ee) + +!*** variables and functions from other modules ***! +use prec, only: pInt, & + pReal +use debug, only: debugger, & + selectiveDebugger, & + verboseDebugger, & + debug_e, & + debug_i, & + debug_g, & + debug_StateLoopDistribution +use FEsolving, only: FEsolving_execElem, & + FEsolving_execIP +use mesh, only: mesh_element, & + mesh_NcpElems +use material, only: homogenization_Ngrains +use constitutive, only: constitutive_sizeDotState, & + constitutive_state, & + constitutive_subState0, & + constitutive_dotState, & + constitutive_collectDotState, & + constitutive_dotTemperature, & + constitutive_microstructure + +implicit none + +!*** input variables ***! +integer(pInt), intent(in) :: mode ! mode of calculation; 1: central solution, 2: stiffness (by perturbation) +integer(pInt), optional, intent(in):: ee, & ! element index + ii, & ! integration point index + gg ! grain index + +!*** output variables ***! + +!*** local variables ***! +integer(pInt) e, & ! element index in element loop + i, & ! integration point index in ip loop + g, & ! grain index in grain loop + n, & + mySizeDotState +integer(pInt), dimension(2) :: eIter ! bounds for element iteration +integer(pInt), dimension(2,mesh_NcpElems) :: iIter, & ! bounds for ip iteration + gIter ! bounds for grain iteration +logical singleRun ! flag indicating computation for single (g,i,e) triple + + +if (present(ee) .and. present(ii) .and. present(gg)) then + eIter = ee + iIter(:,ee) = ii + gIter(:,ee) = gg + singleRun = .true. +else + eIter = FEsolving_execElem(1:2) + do e = eIter(1),eIter(2) + iIter(:,e) = FEsolving_execIP(1:2,e) + gIter(:,e) = (/1,homogenization_Ngrains(mesh_element(3,e))/) + enddo + singleRun = .false. +endif + + +! --- UPDATE DEPENDENT STATES --- + +!$OMP PARALLEL DO + do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, & + crystallite_Fp, crystallite_disorientation(:,:,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states + endif + enddo; enddo; enddo +!$OMPEND PARALLEL DO + + +! --- DOT STATE AND TEMPERATURE --- + +!$OMP PARALLEL DO + do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + selectiveDebugger = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g .and. mode == 1) + call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & + crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & + crystallite_disorientation(:,:,g,i,e), crystallite_subdt(g,i,e), g,i,e) + crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), & + crystallite_Temperature(g,i,e),g,i,e) + if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState + .or. crystallite_dotTemperature(g,i,e)/=crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature + if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local... + !$OMP CRITICAL + crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped + !$OMPEND CRITICAL + else ! if broken local... + crystallite_todo(g,i,e) = .false. ! ... skip this one next time + endif + endif + endif + enddo; enddo; enddo +!$OMPEND PARALLEL DO + + +! --- UPDATE STATE AND TEMPERATURE --- + +!$OMP PARALLEL DO + do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode == 1) + mySizeDotState = constitutive_sizeDotState(g,i,e) + constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_subState0(g,i,e)%p(1:mySizeDotState) & + + constitutive_dotState(g,i,e)%p(1:mySizeDotState) * crystallite_subdt(g,i,e) + crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) & + + crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e) + + if (verboseDebugger .and. selectiveDebugger) then + !$OMP CRITICAL (write2out) + write(6,*) '::: updateState',g,i,e + write(6,*) + write(6,'(a,/,12(e14.8,x))') 'updateState: dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState) + write(6,*) + write(6,'(a,/,12(e14.8,x))') 'updateState: new state', constitutive_state(g,i,e)%p(1:mySizeDotState) + write(6,*) + !$OMPEND CRITICAL (write2out) + endif + endif + enddo; enddo; enddo +!$OMPEND PARALLEL DO + + +! --- UPDATE DEPENDENT STATES --- + +!$OMP PARALLEL DO + do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, & + crystallite_Fp, crystallite_disorientation(:,:,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states + endif + enddo; enddo; enddo +!$OMPEND PARALLEL DO + + +! --- STRESS INTEGRATION --- + +!$OMP PARALLEL DO + do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + selectiveDebugger = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g .and. mode == 1) + if (crystallite_todo(g,i,e)) then + if (crystallite_integrateStress(mode,g,i,e)) then + crystallite_converged(g,i,e) = .true. + !$OMP CRITICAL (distributionState) + debug_StateLoopDistribution(1,mode) = debug_StateLoopDistribution(1,mode) + 1 + !$OMPEND CRITICAL (distributionState) + else ! broken stress integration + if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local... + !$OMP CRITICAL + crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped + !$OMPEND CRITICAL + endif + endif + endif + enddo; enddo; enddo +!$OMPEND PARALLEL DO + + +! --- CHECK NON-LOCAL CONVERGENCE --- + +crystallite_todo = .false. ! done with integration +if ( .not. (mode == 2 .and. singleRun) & ! except for local stiffness calculation: + .and. any(.not. crystallite_converged .and. .not. crystallite_localConstitution)) then ! any non-local not yet converged (or broken)... + crystallite_converged = crystallite_converged .and. crystallite_localConstitution ! ...restart all non-local as not converged +endif + +endsubroutine + + + +!******************************************************************** +! integrate stress, state and Temperature with +! adaptive 1st order explicit Euler method +! using Fixed Point Iteration to adapt the stepsize +!******************************************************************** +subroutine crystallite_integrateStateFPI(mode,gg,ii,ee) + +!*** variables and functions from other modules ***! +use prec, only: pInt, & + pReal +use debug, only: debugger, & + selectiveDebugger, & + verboseDebugger, & + debug_e, & + debug_i, & + debug_g, & + debug_StateLoopDistribution +use numerics, only: nState +use FEsolving, only: FEsolving_execElem, & + FEsolving_execIP +use mesh, only: mesh_element, & + mesh_NcpElems +use material, only: homogenization_Ngrains +use constitutive, only: constitutive_sizeDotState, & + constitutive_state, & + constitutive_dotState, & + constitutive_collectDotState, & + constitutive_dotTemperature, & + constitutive_microstructure, & + constitutive_previousDotState, & + constitutive_previousDotState2 + +implicit none + +!*** input variables ***! +integer(pInt), intent(in) :: mode ! mode of calculation; 1: central solution, 2: stiffness (by perturbation) +integer(pInt), optional, intent(in):: ee, & ! element index + ii, & ! integration point index + gg ! grain index + +!*** output variables ***! + +!*** local variables ***! +integer(pInt) NiterationState, & ! number of iterations in state loop + e, & ! element index in element loop + i, & ! integration point index in ip loop + g ! grain index in grain loop +integer(pInt), dimension(2) :: eIter ! bounds for element iteration +integer(pInt), dimension(2,mesh_NcpElems) :: iIter, & ! bounds for ip iteration + gIter ! bounds for grain iteration +real(pReal) dot_prod12, & + dot_prod22 +logical singleRun ! flag indicating computation for single (g,i,e) triple + + +if (present(ee) .and. present(ii) .and. present(gg)) then + eIter = ee + iIter(:,ee) = ii + gIter(:,ee) = gg + singleRun = .true. +else + eIter = FEsolving_execElem(1:2) + do e = eIter(1),eIter(2) + iIter(:,e) = FEsolving_execIP(1:2,e) + gIter(:,e) = (/1,homogenization_Ngrains(mesh_element(3,e))/) + enddo + singleRun = .false. +endif + + +! --+>> PREGUESS FOR STATE <<+-- + +! --- UPDATE DEPENDENT STATES --- + +!$OMP PARALLEL DO + do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, & + crystallite_Fp, crystallite_disorientation(:,:,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states + endif + enddo; enddo; enddo +!$OMPEND PARALLEL DO + + +! --- DOT STATES --- + +!$OMP PARALLEL DO + do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + selectiveDebugger = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g .and. mode == 1) + if (crystallite_todo(g,i,e)) then + constitutive_previousDotState2(g,i,e)%p = 0.0_pReal + constitutive_previousDotState(g,i,e)%p = 0.0_pReal + call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & + crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & + crystallite_disorientation(:,:,g,i,e), crystallite_subdt(g,i,e), g, i, e) + endif + enddo; enddo; enddo +!$OMPEND PARALLEL DO + + +! --- STATE & TEMPERATURE UPDATE --- + +crystallite_statedamper = 1.0_pReal +!$OMP PARALLEL DO + do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode == 1) + if (crystallite_todo(g,i,e)) then + crystallite_stateConverged(g,i,e) = crystallite_updateState(g,i,e) ! update state + crystallite_temperatureConverged(g,i,e) = crystallite_updateTemperature(g,i,e) ! update temperature + if ( .not. crystallite_localConstitution(g,i,e) .and. .not. crystallite_todo(g,i,e) ) then ! if updateState or updateTemperature signals broken non-local... + !$OMP CRITICAL + crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped + !$OMPEND CRITICAL + endif + endif + enddo; enddo; enddo +!$OMPEND PARALLEL DO + + +! --+>> STATE LOOP <<+-- + +NiterationState = 0_pInt + +do while (any(crystallite_todo) .and. NiterationState < nState ) ! convergence loop for crystallite + + NiterationState = NiterationState + 1_pInt + + + ! --- STRESS INTEGRATION --- + + !$OMP PARALLEL DO + do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + selectiveDebugger = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g .and. mode == 1) + if (crystallite_todo(g,i,e)) then + crystallite_todo(g,i,e) = crystallite_integrateStress(mode,g,i,e) + if ( .not. crystallite_localConstitution(g,i,e) .and. .not. crystallite_todo(g,i,e)) then ! if broken non-local... + !$OMP CRITICAL + crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped + !$OMPEND CRITICAL + endif + endif + enddo; enddo; enddo + !$OMPEND PARALLEL DO + + if (verboseDebugger .and. mode == 1) then + !$OMP CRITICAL (write2out) + write(6,*) count(crystallite_todo(:,:,:)),'grains todo after stress integration' + !$OMPEND CRITICAL (write2out) + endif + + + ! --- DOT STATES --- + + !$OMP PARALLEL DO + do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + selectiveDebugger = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g .and. mode == 1) + if (crystallite_todo(g,i,e)) then + constitutive_previousDotState2(g,i,e)%p = constitutive_previousDotState(g,i,e)%p ! wind forward dotStates + constitutive_previousDotState(g,i,e)%p = constitutive_dotState(g,i,e)%p + call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & + crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & + crystallite_disorientation(:,:,g,i,e), crystallite_subdt(g,i,e), g, i, e) + endif + enddo; enddo; enddo + !$OMPEND PARALLEL DO + + + ! --- STATE & TEMPERATURE UPDATE --- + + crystallite_statedamper = 1.0_pReal + !$OMP PARALLEL DO + do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode == 1) + if (crystallite_todo(g,i,e)) then + + ! --- state damper --- + + dot_prod12 = dot_product( constitutive_dotState(g,i,e)%p - constitutive_previousDotState(g,i,e)%p, & + constitutive_previousDotState(g,i,e)%p - constitutive_previousDotState2(g,i,e)%p ) + dot_prod22 = dot_product( constitutive_previousDotState(g,i,e)%p - constitutive_previousDotState2(g,i,e)%p, & + constitutive_previousDotState(g,i,e)%p - constitutive_previousDotState2(g,i,e)%p ) + if ( dot_prod22 > 0.0_pReal & + .and. ( dot_prod12 < 0.0_pReal & + .or. dot_product(constitutive_dotState(g,i,e)%p, constitutive_previousDotState(g,i,e)%p) < 0.0_pReal) ) & + crystallite_statedamper(g,i,e) = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22) + + ! --- updates --- + + crystallite_stateConverged(g,i,e) = crystallite_updateState(g,i,e) ! update state + crystallite_temperatureConverged(g,i,e) = crystallite_updateTemperature(g,i,e) ! update temperature + crystallite_converged(g,i,e) = crystallite_stateConverged(g,i,e) .and. crystallite_temperatureConverged(g,i,e) + if ( .not. crystallite_localConstitution(g,i,e) .and. .not. crystallite_todo(g,i,e)) then ! if updateState or updateTemperature signals broken non-local... + !$OMP CRITICAL + crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped + !$OMPEND CRITICAL + elseif (crystallite_converged(g,i,e)) then + !$OMP CRITICAL (distributionState) + debug_StateLoopDistribution(NiterationState,mode) = debug_StateLoopDistribution(NiterationState,mode) + 1 + !$OMPEND CRITICAL (distributionState) + endif + endif + enddo; enddo; enddo + !$OMPEND PARALLEL DO + + + ! --- UPDATE DEPENDENT STATES --- + + !$OMP PARALLEL DO + do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, & + crystallite_Fp, crystallite_disorientation(:,:,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states + endif + enddo; enddo; enddo + !$OMPEND PARALLEL DO + + if (verboseDebugger .and. mode == 1) then + !$OMP CRITICAL (write2out) + write(6,*) count(crystallite_converged(:,:,:)),'grains converged after state integration no.', NiterationState + write(6,*) + !$OMPEND CRITICAL (write2out) + endif + + + ! --- CONVERGENCE CHECK --- + + if ( .not. (mode == 2 .and. singleRun) & ! except for local stiffness calculation: + .and. any(.not. crystallite_converged .and. .not. crystallite_localConstitution)) then ! any non-local not yet converged (or broken)... + crystallite_converged = crystallite_converged .and. crystallite_localConstitution ! ...restart all non-local as not converged + endif + + crystallite_todo = crystallite_todo .and. .not. crystallite_converged ! skip all converged + + if (verboseDebugger .and. mode == 1) then + !$OMP CRITICAL (write2out) + write(6,*) count(crystallite_converged(:,:,:)),'grains converged after non-local check' + write(6,*) count(crystallite_todo(:,:,:)),'grains todo after state integration no.', NiterationState + write(6,*) + !$OMPEND CRITICAL (write2out) + endif + +enddo ! crystallite convergence loop + +endsubroutine + + + !******************************************************************** ! update the internal state of the constitutive law ! and tell whether state has converged !******************************************************************** - function crystallite_updateState(& - g,& ! grain number - i,& ! integration point number - e & ! element number - ) - - !*** variables and functions from other modules ***! - use prec, only: pReal, & - pInt, & - pLongInt - use numerics, only: rTol_crystalliteState - use constitutive, only: constitutive_dotState, & - constitutive_previousDotState, & - constitutive_sizeDotState, & - constitutive_subState0, & - constitutive_state, & - constitutive_relevantState, & - constitutive_microstructure - use debug, only: debugger, & - selectiveDebugger, & - verboseDebugger - use FEsolving, only: cycleCounter, theInc - - !*** input variables ***! - integer(pInt), intent(in):: e, & ! element index - i, & ! integration point index - g ! grain index - - !*** output variables ***! - logical crystallite_updateState ! flag indicating if integration suceeded +function crystallite_updateState(g,i,e) - !*** local variables ***! - real(pReal), dimension(constitutive_sizeDotState(g,i,e)) :: residuum ! residuum from evolution of microstructure - integer(pInt) mySize +!*** variables and functions from other modules ***! +use prec, only: pReal, & + pInt, & + pLongInt +use numerics, only: rTol_crystalliteState +use constitutive, only: constitutive_dotState, & + constitutive_previousDotState, & + constitutive_sizeDotState, & + constitutive_subState0, & + constitutive_state, & + constitutive_relevantState, & + constitutive_microstructure +use debug, only: debugger, & + selectiveDebugger, & + verboseDebugger - - mySize = constitutive_sizeDotState(g,i,e) +!*** input variables ***! +integer(pInt), intent(in):: e, & ! element index + i, & ! integration point index + g ! grain index - ! correct my dotState - constitutive_dotState(g,i,e)%p(1:mySize) = constitutive_dotState(g,i,e)%p(1:mySize) * crystallite_statedamper(g,i,e) & - + constitutive_previousDotState(g,i,e)%p(1:mySize) * & - (1.0_pReal-crystallite_statedamper(g,i,e)) - ! calculate the residuum - residuum = constitutive_state(g,i,e)%p(1:mySize) - constitutive_subState0(g,i,e)%p(1:mySize) & - - constitutive_dotState(g,i,e)%p(1:mySize) * crystallite_subdt(g,i,e) - - if (any(residuum/=residuum)) then ! if NaN occured then return without changing the state... - crystallite_updateState = .false. ! ...indicate state update failed - crystallite_todo(g,i,e) = .false. ! ...no need to calculate any further - if (verboseDebugger) then - !$OMP CRITICAL (write2out) - write(6,*) '::: updateState encountered NaN',g,i,e - !$OMPEND CRITICAL (write2out) - endif - return - endif - - ! update the microstructure - constitutive_state(g,i,e)%p(1:mySize) = constitutive_state(g,i,e)%p(1:mySize) - residuum - call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, crystallite_Fp, & - crystallite_disorientation(:,:,g,i,e), g, i, e) - - - ! setting flag to true if state is below relative tolerance, otherwise set it to false - crystallite_updateState = all( constitutive_state(g,i,e)%p(1:mySize) < constitutive_relevantState(g,i,e)%p(1:mySize) & +!*** output variables ***! +logical crystallite_updateState ! flag indicating if integration suceeded + +!*** local variables ***! +real(pReal), dimension(constitutive_sizeDotState(g,i,e)) :: residuum ! residuum from evolution of microstructure +integer(pInt) mySize + + +mySize = constitutive_sizeDotState(g,i,e) + +! correct my dotState +constitutive_dotState(g,i,e)%p(1:mySize) = constitutive_dotState(g,i,e)%p(1:mySize) * crystallite_statedamper(g,i,e) & + + constitutive_previousDotState(g,i,e)%p(1:mySize) * (1.0_pReal-crystallite_statedamper(g,i,e)) + +residuum = constitutive_state(g,i,e)%p(1:mySize) - constitutive_subState0(g,i,e)%p(1:mySize) & + - constitutive_dotState(g,i,e)%p(1:mySize) * crystallite_subdt(g,i,e) + +if (any(residuum/=residuum)) then ! if NaN occured then return without changing the state... + crystallite_updateState = .false. ! ...indicate state update failed + crystallite_todo(g,i,e) = .false. ! ...no need to calculate any further + if (verboseDebugger) then + !$OMP CRITICAL (write2out) + write(6,*) '::: updateState encountered NaN',g,i,e + !$OMPEND CRITICAL (write2out) + endif + return +endif + +constitutive_state(g,i,e)%p(1:mySize) = constitutive_state(g,i,e)%p(1:mySize) - residuum + +! setting flag to true if state is below relative tolerance, otherwise set it to false +crystallite_updateState = all( constitutive_state(g,i,e)%p(1:mySize) < constitutive_relevantState(g,i,e)%p(1:mySize) & .or. abs(residuum) < rTol_crystalliteState*abs(constitutive_state(g,i,e)%p(1:mySize)) ) - if (verboseDebugger .and. selectiveDebugger) then - !$OMP CRITICAL (write2out) - if (crystallite_updateState) then - write(6,*) '::: updateState converged',g,i,e - else - write(6,*) '::: updateState did not converge',g,i,e - endif - write(6,*) - write(6,'(a,f6.1)') 'updateState: crystallite_statedamper',crystallite_statedamper(g,i,e) - write(6,*) - write(6,'(a,/,12(e14.8,x))') 'updateState: dotState',constitutive_dotState(g,i,e)%p(1:mySize) - write(6,*) - write(6,'(a,/,12(e14.8,x))') 'updateState: new state',constitutive_state(g,i,e)%p(1:mySize) - write(6,*) - write(6,'(a,/,12(f12.1,x))') 'updateState: resid tolerance',abs(residuum/rTol_crystalliteState/ & - constitutive_state(g,i,e)%p(1:mySize)) - write(6,*) - !$OMPEND CRITICAL (write2out) - endif - return - endfunction +if (verboseDebugger .and. selectiveDebugger) then + !$OMP CRITICAL (write2out) + if (crystallite_updateState) then + write(6,*) '::: updateState converged',g,i,e + else + write(6,*) '::: updateState did not converge',g,i,e + endif + write(6,*) + write(6,'(a,f6.1)') 'updateState: crystallite_statedamper',crystallite_statedamper(g,i,e) + write(6,*) + write(6,'(a,/,12(e14.8,x))') 'updateState: dotState',constitutive_dotState(g,i,e)%p(1:mySize) + write(6,*) + write(6,'(a,/,12(e14.8,x))') 'updateState: new state',constitutive_state(g,i,e)%p(1:mySize) + write(6,*) + write(6,'(a,/,12(f12.1,x))') 'updateState: resid tolerance',abs(residuum / rTol_crystalliteState & + / constitutive_state(g,i,e)%p(1:mySize)) + write(6,*) + !$OMPEND CRITICAL (write2out) +endif + +endfunction + !******************************************************************** @@ -1180,9 +2255,7 @@ endsubroutine pLongInt use numerics, only: rTol_crystalliteTemperature use constitutive, only: constitutive_dotTemperature - use debug, only: debugger, & - debug_cumDotTemperatureCalls, & - debug_cumDotTemperatureTicks + use debug, only: debugger !*** input variables ***! integer(pInt), intent(in):: e, & ! element index @@ -1194,24 +2267,16 @@ endsubroutine !*** local variables ***! real(pReal) residuum ! residuum from evolution of temperature - integer(pLongInt) tick, & - tock, & - tickrate, & - maxticks ! calculate the residuum - call system_clock(count=tick,count_rate=tickrate,count_max=maxticks) residuum = crystallite_Temperature(g,i,e) - crystallite_subTemperature0(g,i,e) - & crystallite_subdt(g,i,e) * & constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e),crystallite_Temperature(g,i,e),g,i,e) - call system_clock(count=tock,count_rate=tickrate,count_max=maxticks) - debug_cumDotTemperatureCalls = debug_cumDotTemperatureCalls + 1_pInt - debug_cumDotTemperatureTicks = debug_cumDotTemperatureTicks + tock-tick - if (tock < tick) debug_cumDotTemperatureTicks = debug_cumDotTemperatureTicks + maxticks ! if NaN occured then return without changing the state if (residuum/=residuum) then crystallite_updateTemperature = .false. ! indicate update failed + crystallite_todo(g,i,e) = .false. ! ...no need to calculate any further !$OMP CRITICAL (write2out) write(6,*) '::: updateTemperature encountered NaN',g,i,e !$OMPEND CRITICAL (write2out) @@ -1240,7 +2305,9 @@ endsubroutine mode, & ! 1: central solution, 2: stiffness (by perturbation) g,& ! grain number i,& ! integration point number - e) ! element number + e,& ! element number + fraction & + ) !*** variables and functions from other modules ***! @@ -1281,7 +2348,8 @@ endsubroutine e, & ! element index i, & ! integration point index g ! grain index - + real(pReal), optional, intent(in) :: fraction ! fraction of timestep + !*** output variables ***! logical crystallite_integrateStress ! flag indicating if integration suceeded @@ -1312,7 +2380,8 @@ endsubroutine real(pReal) p_hydro, & ! volumetric part of 2nd Piola-Kirchhoff Stress det, & ! determinant leapfrog, & ! acceleration factor for Newton-Raphson scheme - maxleap ! maximum acceleration factor + maxleap, & ! maximum acceleration factor + dt ! time increment logical error ! flag indicating an error integer(pInt) NiterationStress, & ! number of stress integrations dummy, & @@ -1331,8 +2400,16 @@ endsubroutine ! be pessimistic crystallite_integrateStress = .false. + ! only integrate over fraction of timestep? + if (present(fraction)) then + dt = crystallite_subdt(g,i,e) * fraction + Fg_new = crystallite_subF0(:,:,g,i,e) + (crystallite_subF(:,:,g,i,e) - crystallite_subF0(:,:,g,i,e)) * fraction + else + dt = crystallite_subdt(g,i,e) + Fg_new = crystallite_subF(:,:,g,i,e) + endif + ! feed local variables - Fg_new = crystallite_subF(:,:,g,i,e) Fp_current = crystallite_subFp0(:,:,g,i,e) Tstar_v = crystallite_Tstar_v(:,g,i,e) Lpguess_old = crystallite_Lp(:,:,g,i,e) ! consider present Lp good (i.e. worth remembering) ... @@ -1381,7 +2458,7 @@ LpLoop: do return endif - B = math_I3 - crystallite_subdt(g,i,e)*Lpguess + B = math_I3 - dt*Lpguess BT = transpose(B) AB = math_mul33x33(A,B) BTA = math_mul33x33(BT,A) @@ -1413,8 +2490,8 @@ LpLoop: do ! Check for convergence of loop if (.not.(any(residuum/=residuum)) .and. & ! exclude any NaN in residuum ( maxval(abs(residuum)) < aTol_crystalliteStress .or. & ! below absolute tolerance .or. - ( any(abs(crystallite_subdt(g,i,e)*Lpguess) > relevantStrain) .and. & ! worth checking? .and. - maxval(abs(residuum/Lpguess), abs(crystallite_subdt(g,i,e)*Lpguess) > relevantStrain) < rTol_crystalliteStress & ! below relative tolerance + ( any(abs(dt*Lpguess) > relevantStrain) .and. & ! worth checking? .and. + maxval(abs(residuum/Lpguess), abs(dt*Lpguess) > relevantStrain) < rTol_crystalliteStress & ! below relative tolerance ) & ) & ) & @@ -1463,7 +2540,7 @@ LpLoop: do ! forall (h=1:3,j=1:3,k=1:3,l=1:3,m=1:3) & dTdLp(3*(h-1)+j,3*(k-1)+l) = dTdLp(3*(h-1)+j,3*(k-1)+l) + C(h,j,l,m)*AB(k,m)+C(h,j,m,l)*BTA(m,k) enddo; enddo; enddo; enddo; enddo - dTdLp = -0.5_pReal*crystallite_subdt(g,i,e)*dTdLp + dTdLp = -0.5_pReal*dt*dTdLp dRdLp = math_identity2nd(9) - math_mul99x99(dLpdT_constitutive,dTdLp) invdRdLp = 0.0_pReal call math_invert(9,dRdLp,invdRdLp,dummy,error) ! invert dR/dLp --> dLp/dR @@ -1738,16 +2815,16 @@ function crystallite_postResults(& c = c + 1_pInt case ('orientation') crystallite_postResults(c+1:c+4) = & - crystallite_orientation(:,g,i,e) ! grain orientation as quaternion + crystallite_orientation(:,g,i,e) ! grain orientation as quaternion c = c + 4_pInt case ('eulerangles') crystallite_postResults(c+1:c+3) = inDeg * & - math_QuaternionToEuler(crystallite_orientation(:,g,i,e)) ! grain orientation as Euler angles in degree + math_QuaternionToEuler(crystallite_orientation(:,g,i,e)) ! grain orientation as Euler angles in degree c = c + 3_pInt case ('grainrotation') crystallite_postResults(c+1:c+4) = & - math_QuaternionToAxisAngle(crystallite_rotation(1:4,g,i,e)) ! grain rotation away from initial orientation as axis-angle - crystallite_postResults(c+4) = inDeg * crystallite_postResults(c+4) ! angle in degree + math_QuaternionToAxisAngle(crystallite_rotation(1:4,g,i,e)) ! grain rotation away from initial orientation as axis-angle + crystallite_postResults(c+4) = inDeg * crystallite_postResults(c+4) ! angle in degree c = c + 4_pInt case ('defgrad','f') mySize = 9_pInt @@ -1779,7 +2856,7 @@ function crystallite_postResults(& crystallite_postResults(c+1) = constitutive_sizePostResults(g,i,e); c = c+1_pInt ! size of constitutive results crystallite_postResults(c+1:c+constitutive_sizePostResults(g,i,e)) = & - constitutive_postResults(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, crystallite_Fp, & + constitutive_postResults(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, crystallite_Fp, & crystallite_Temperature(g,i,e), crystallite_disorientation(:,:,g,i,e), dt, & crystallite_subdt(g,i,e), g, i, e) c = c + constitutive_sizePostResults(g,i,e) diff --git a/code/debug.config b/code/debug.config index 46ea3c8c4..c6414c580 100644 --- a/code/debug.config +++ b/code/debug.config @@ -3,7 +3,7 @@ debug 1 # >0 true to switch on debugging verbose 1 # >0 true to switch on verbose output -selective 1 # >0 true to switch on e,i,g seelctive debugging -element 123 # selected element for debugging (synonymous: "el", "e") -ip 3 # selected integration point for debugging (synonymous: "integrationpoint", "i") -grain 5 # selected grain at ip for debugging (synonymous: "gr", "g") +selective 1 # >0 true to switch on e,i,g selective debugging +element 1 # selected element for debugging (synonymous: "el", "e") +ip 1 # selected integration point for debugging (synonymous: "integrationpoint", "i") +grain 1 # selected grain at ip for debugging (synonymous: "gr", "g") diff --git a/code/debug.f90 b/code/debug.f90 index baf104261..882dbdc29 100644 --- a/code/debug.f90 +++ b/code/debug.f90 @@ -9,8 +9,7 @@ integer(pInt), dimension(:,:), allocatable :: debug_StressLoopDistribution integer(pInt), dimension(:,:), allocatable :: debug_LeapfrogBreakDistribution - integer(pInt), dimension(:), allocatable :: debug_CrystalliteStateLoopDistribution - integer(pInt), dimension(:), allocatable :: debug_StiffnessStateLoopDistribution + integer(pInt), dimension(:,:), allocatable :: debug_StateLoopDistribution integer(pInt), dimension(:), allocatable :: debug_CrystalliteLoopDistribution integer(pInt), dimension(:), allocatable :: debug_MaterialpointStateLoopDistribution integer(pInt), dimension(:), allocatable :: debug_MaterialpointLoopDistribution @@ -70,8 +69,7 @@ subroutine debug_init() allocate(debug_StressLoopDistribution(nStress,2)) ; debug_StressLoopDistribution = 0_pInt allocate(debug_LeapfrogBreakDistribution(nStress,2)) ; debug_LeapfrogBreakDistribution = 0_pInt - allocate(debug_CrystalliteStateLoopDistribution(nState)) ; debug_CrystalliteStateLoopDistribution = 0_pInt - allocate(debug_StiffnessStateLoopDistribution(nState)) ; debug_StiffnessStateLoopDistribution = 0_pInt + allocate(debug_StateLoopDistribution(nState,2)) ; debug_StateLoopDistribution = 0_pInt allocate(debug_CrystalliteLoopDistribution(nCryst+1)) ; debug_CrystalliteLoopDistribution = 0_pInt allocate(debug_MaterialpointStateLoopDistribution(nMPstate)) ; debug_MaterialpointStateLoopDistribution = 0_pInt allocate(debug_MaterialpointLoopDistribution(nHomog+1)) ; debug_MaterialpointLoopDistribution = 0_pInt @@ -141,8 +139,7 @@ subroutine debug_reset() debug_StressLoopDistribution = 0_pInt ! initialize debugging data debug_LeapfrogBreakDistribution = 0_pInt - debug_CrystalliteStateLoopDistribution = 0_pInt - debug_StiffnessStateLoopDistribution = 0_pInt + debug_StateLoopDistribution = 0_pInt debug_CrystalliteLoopDistribution = 0_pInt debug_MaterialpointStateLoopDistribution = 0_pInt debug_MaterialpointLoopDistribution = 0_pInt @@ -216,12 +213,14 @@ endsubroutine write(6,*) write(6,*) 'distribution_CrystalliteStateLoop :' do i=1,nState - if (debug_CrystalliteStateLoopDistribution(i) /= 0) then - integral = integral + i*debug_CrystalliteStateLoopDistribution(i) - write(6,'(i25,x,i10)') i,debug_CrystalliteStateLoopDistribution(i) + if (any(debug_StateLoopDistribution(i,:) /= 0)) then + integral = integral + i*debug_StateLoopDistribution(i,1) + i*debug_StateLoopDistribution(i,2) + write(6,'(i25,x,i10,12x,i10)') i,debug_StateLoopDistribution(i,1),debug_StateLoopDistribution(i,2) endif enddo - write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_CrystalliteStateLoopDistribution) + write(6,'(a15,i10,x,i10,12x,i10)') ' total',integral,& + sum(debug_StateLoopDistribution(:,1)), & + sum(debug_StateLoopDistribution(:,2)) integral = 0_pInt write(6,*) @@ -237,19 +236,7 @@ endsubroutine endif enddo write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_CrystalliteLoopDistribution) - - integral = 0_pInt - write(6,*) - write(6,*) 'distribution_StiffnessStateLoop :' - do i=1,nState - if (debug_StiffnessStateLoopDistribution(i) /= 0) then - integral = integral + i*debug_StiffnessStateLoopDistribution(i) - write(6,'(i25,x,i10)') i,debug_StiffnessStateLoopDistribution(i) - endif - enddo - write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_StiffnessStateLoopDistribution) -!* Material point loop counter <<>> integral = 0_pInt write(6,*) write(6,*) diff --git a/code/homogenization.f90 b/code/homogenization.f90 index df9bb97e2..1808a9bd5 100644 --- a/code/homogenization.f90 +++ b/code/homogenization.f90 @@ -229,6 +229,7 @@ subroutine materialpoint_stressAndItsTangent(& stepIncreaseHomog, & nHomog, & nMPstate + use math, only: math_det3x3 use FEsolving, only: FEsolving_execElem, & FEsolving_execIP, & terminallyIll @@ -243,6 +244,8 @@ subroutine materialpoint_stressAndItsTangent(& crystallite_Fp, & crystallite_Lp0, & crystallite_Lp, & + crystallite_dPdF, & + crystallite_dPdF0, & crystallite_Tstar0_v, & crystallite_Tstar_v, & crystallite_partionedTemperature0, & @@ -250,6 +253,7 @@ subroutine materialpoint_stressAndItsTangent(& crystallite_partionedF, & crystallite_partionedFp0, & crystallite_partionedLp0, & + crystallite_partioneddPdF0, & crystallite_partionedTstar0_v, & crystallite_dt, & crystallite_requested, & @@ -296,6 +300,7 @@ subroutine materialpoint_stressAndItsTangent(& crystallite_partionedTemperature0(1:myNgrains,i,e) = materialpoint_Temperature(i,e) ! ...temperatures crystallite_partionedFp0(:,:,1:myNgrains,i,e) = crystallite_Fp0(:,:,1:myNgrains,i,e) ! ...plastic def grads crystallite_partionedLp0(:,:,1:myNgrains,i,e) = crystallite_Lp0(:,:,1:myNgrains,i,e) ! ...plastic velocity grads + crystallite_partioneddPdF0(:,:,:,:,1:myNgrains,i,e) = crystallite_dPdF0(:,:,:,:,1:myNgrains,i,e) ! ...stiffness crystallite_partionedF0(:,:,1:myNgrains,i,e) = crystallite_F0(:,:,1:myNgrains,i,e) ! ...def grads crystallite_partionedTstar0_v(:,1:myNgrains,i,e)= crystallite_Tstar0_v(:,1:myNgrains,i,e) ! ...2nd PK stress @@ -348,6 +353,7 @@ subroutine materialpoint_stressAndItsTangent(& crystallite_partionedF0(:,:,1:myNgrains,i,e) = crystallite_partionedF(:,:,1:myNgrains,i,e) ! ...def grads crystallite_partionedFp0(:,:,1:myNgrains,i,e) = crystallite_Fp(:,:,1:myNgrains,i,e) ! ...plastic def grads crystallite_partionedLp0(:,:,1:myNgrains,i,e) = crystallite_Lp(:,:,1:myNgrains,i,e) ! ...plastic velocity grads + crystallite_partioneddPdF0(:,:,:,:,1:myNgrains,i,e) = crystallite_dPdF(:,:,:,:,1:myNgrains,i,e)! ...stiffness crystallite_partionedTstar0_v(:,1:myNgrains,i,e) = crystallite_Tstar_v(:,1:myNgrains,i,e) ! ...2nd PK stress forall (g = 1:myNgrains) constitutive_partionedState0(g,i,e)%p = constitutive_state(g,i,e)%p ! ...microstructures if (homogenization_sizeState(i,e) > 0_pInt) & @@ -382,6 +388,7 @@ subroutine materialpoint_stressAndItsTangent(& ! ...initial def grad unchanged crystallite_Fp(:,:,1:myNgrains,i,e) = crystallite_partionedFp0(:,:,1:myNgrains,i,e) ! ...plastic def grads crystallite_Lp(:,:,1:myNgrains,i,e) = crystallite_partionedLp0(:,:,1:myNgrains,i,e) ! ...plastic velocity grads + crystallite_dPdF(:,:,:,:,1:myNgrains,i,e) = crystallite_partioneddPdF0(:,:,:,:,1:myNgrains,i,e) ! ...stiffness crystallite_Tstar_v(:,1:myNgrains,i,e) = crystallite_partionedTstar0_v(:,1:myNgrains,i,e) ! ...2nd PK stress forall (g = 1:myNgrains) constitutive_state(g,i,e)%p = constitutive_partionedState0(g,i,e)%p ! ...microstructures if (homogenization_sizeState(i,e) > 0_pInt) & @@ -478,13 +485,13 @@ subroutine materialpoint_stressAndItsTangent(& enddo ! homogenization convergence loop - NiterationHomog = NiterationHomog +1_pInt + NiterationHomog = NiterationHomog + 1_pInt enddo ! cutback loop if (.not. terminallyIll ) then - + call crystallite_orientations() ! calculate crystal orientations !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed diff --git a/code/lattice.f90 b/code/lattice.f90 index 3e353e7ef..e446ac914 100644 --- a/code/lattice.f90 +++ b/code/lattice.f90 @@ -26,6 +26,9 @@ integer(pInt), parameter :: lattice_maxNslip = 48 ! max # of sli integer(pInt), parameter :: lattice_maxNtwin = 24 ! max # of twin systems over lattice structures integer(pInt), parameter :: lattice_maxNinteraction = 20 ! max # of interaction types (in hardening matrix part) +integer(pInt), parameter, dimension(3) :: lattice_symmetryTypes =(/1, 1, 2/) ! maps crystal structures to symmetry tpyes + + integer(pInt), pointer, dimension(:,:) :: interactionSlipSlip, & interactionSlipTwin, & interactionTwinSlip, & @@ -429,37 +432,37 @@ integer(pInt), allocatable, dimension(:,:,:) :: lattice_interactionSlipSlip, & integer(pInt), parameter :: lattice_hex_Ntwin = 24 ! sum(lattice_hex_NtwinSystem) integer(pInt) :: lattice_hex_Nstructure = 0_pInt - !* sorted by A. Alankar & P. Eisenlohr + !* sorted by YJ.Ro and Philip real(pReal), dimension(4+4,lattice_hex_Nslip), parameter :: lattice_hex_systemSlip = & reshape((/& ! Basal systems <1120>{0001} (independent of c/a-ratio, Bravais notation (4 coordinate base)) - 2, -1, -1, 0, 0, 0, 0, 1, & !A1 - -1, 2, -1, 0, 0, 0, 0, 1, & !A2 - -1, -1, 2, 0, 0, 0, 0, 1, & !A3 + 2, -1, -1, 0, 0, 0, 0, 1, & + 1, 1, -2, 0, 0, 0, 0, 1, & + -1, 2, -1, 0, 0, 0, 0, 1, & ! 1st type prismatic systems <1120>{1010} (independent of c/a-ratio) - 2, -1, -1, 0, 0, 1, -1, 0, & !B1 - -1, 2, -1, 0, -1, 0, 1, 0, & !C2 - -1, -1, 2, 0, 1, -1, 0, 0, & !D3 -! 1st type 1st order pyramidal systems <1120>{1011} -- plane normals depend on the c/a-ratio - 2, -1, -1, 0, 0, 1, -1, 1, & !E1 - 1, 1, -2, 0, -1, 1, 0, 1, & !F(-3) - -1, 2, -1, 0, -1, 0, 1, 1, & !G2 - -2, 1, 1, 0, 0, -1, 1, 1, & !H(-1) - -1, -1, 2, 0, 1, -1, 0, 1, & !I3 - 1, -2, 1, 0, 1, 0, -1, 1, & !J(-2) + 2, -1, -1, 0, 0, -1, 1, 0, & + 1, 1, 2, 0, 1, -1, 0, 0, & + -1, 2, -1, 0, 1, 0, -1, 0, & +! 1st type 1st order pyramidal systems <1120>{1011} + 2, -1, -1, 0, 0, -1, 1, 1, & + 1, 1, -2, 0, 1, -1, 0, 1, & + -1, 2, -1, 0, 1, 0, -1, 1, & + -2, 1, 1, 0, 0, 1, -1, 1, & + -1, -1, 2, 0, -1, 1, 0, 1, & + 1, -2, 1, 0, -1, 0, 1, 1, & ! pyramidal system: c+a slip <2113>{1011} -- plane normals depend on the c/a-ratio - -1, 2, -1, 3, 0, 1, -1, 1, & !E4 - 1, 1, -2, 3, 0, 1, -1, 1, & !E5 - -2, 1, 1, 3, -1, 1, 0, 1, & !F6 - -1, 2, -1, 3, -1, 1, 0, 1, & !F4 - -1, -1, 2, 3, -1, 0, 1, 1, & !G7 - -2, 1, 1, 3, -1, 0, 1, 1, & !G6 - 1, -2, 1, 3, 0, -1, 1, 1, & !H8 - -1, -1, 2, 3, 0, -1, 1, 1, & !H7 - 2, -1, -1, 3, 1, -1, 0, 1, & !I9 - 1, -2, 1, 3, 1, -1, 0, 1, & !I8 - 1, 1, -2, 3, 1, 0, -1, 1, & !J5 - 2, -1, -1, 3, 1, 0, -1, 1 & !J9 + 2, -1, -1, -3, 1, -1, 0, 1, & + 1, 1, -2, -3, 1, 0, -1, 1, & + -1, 2, -1, -3, 1, -1, 0, 1, & + -2, 1, 1, -3, 1, 0, -1, 1, & + -1, -1, 2, -3, 0, 1, -1, 1, & + 1, -2, 1, -3, 0, -1, 1, 1, & + -2, 1, 1, -3, -1, 0, 1, 1, & + -1, -1, 2, -3, 0, -1, 1, 1, & + 1, -2, 1, -3, 0, 1, -1, 1, & + 2, -1, -1, -3, -1, 1, 0, 1, & + 1, 1, -2, -3, -1, 0, 1, 1, & + -1, 2, -1, -3, -1, 1, 0, 1 & /),(/4+4,lattice_hex_Nslip/)) real(pReal), dimension(4+4,lattice_hex_Ntwin), parameter :: lattice_hex_systemTwin = & @@ -658,30 +661,6 @@ CONTAINS !* - lattice_initializeStructure !**************************************** -pure function lattice_symmetryType(structID) -!************************************** -!* maps structure to symmetry type * -!* fcc(1) and bcc(2) are cubic(1) * -!* hex(3+) is hexagonal(2) * -!************************************** - implicit none - - integer(pInt), intent(in) :: structID - integer(pInt) lattice_symmetryType - - select case(structID) - case (1,2) - lattice_symmetryType = 1_pInt - case (3:) - lattice_symmetryType = 2_pInt - case default - lattice_symmetryType = 0_pInt - end select - - return - -end function - subroutine lattice_init() !************************************** diff --git a/code/material.config b/code/material.config index 9e72b7ca9..ac43c52c2 100644 --- a/code/material.config +++ b/code/material.config @@ -230,9 +230,9 @@ rhoSglScrewPos0 1e11 0 0 0 # Initial positive screw single rhoSglScrewNeg0 1e11 0 0 0 # Initial negative screw single dislocation density in m/m**3 rhoDipEdge0 1e8 0 0 0 # Initial edge dipole dislocation density in m/m**3 rhoDipScrew0 1e8 0 0 0 # Initial screw dipole dislocation density in m/m**3 -r 1e-5 # cutoff radius for dislocation stress in m -v0 1e-4 0 0 0 # prefactor for dislocation velocity -Q0 3e-19 # activation energy for dislocation glide +r 10e-6 # cutoff radius for dislocation stress in m +v0 3500 0 0 0 # maximum dislocation velocity (velocity of sound) +Q0 5e-19 # activation energy for dislocation glide dDipMinEdge 1e-9 0 0 0 # minimum distance for stable edge dipoles in m dDipMinScrew 1e-9 0 0 0 # minimum distance for stable screw dipoles in m lambda0 100 0 0 0 # prefactor for mean free path diff --git a/code/numerics.config b/code/numerics.config index bec167baf..40927970c 100644 --- a/code/numerics.config +++ b/code/numerics.config @@ -7,6 +7,8 @@ iJacoStiffness 1 # frequency of stiffness update iJacoLpresiduum 1 # frequency of Jacobian update of residuum in Lp pert_Fg 1.0e-7 # deformation gradient perturbation for grain tangent pert_method 1 # perturbation method (1 = forward, 2 = backward or 3 = central) +integrator 1 # integration method (1 = Fixed Point Iteration, 2 = Euler, 3 = Adaptive Euler, 4 = classical 4th order Runge-Kutta, 5 = 5th order Runge-Kutta Cash-Karp) +integratorStiffness 1 # integration method used for stiffness (1 = Fixed Point Iteration, 2 = Euler, 3 = Adaptive Euler, 4 = classical 4th order Runge-Kutta, 5 = 5th order Runge-Kutta Cash-Karp) ## crystallite numerical parameters ## nCryst 20 # crystallite loop limit (only for debugging info, loop limit is determined by "subStepMinCryst") diff --git a/code/numerics.f90 b/code/numerics.f90 index cd701b497..51d6dc451 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -14,7 +14,9 @@ integer(pInt) iJacoStiffness, & ! freque nCryst, & ! crystallite loop limit (only for debugging info, loop limit is determined by "subStepMinCryst") nState, & ! state loop limit nStress, & ! stress loop limit - pert_method ! method used in perturbation technique for tangent + pert_method, & ! method used in perturbation technique for tangent + integrator, & ! method used for state integration + integratorStiffness ! method used for stiffness state integration real(pReal) relevantStrain, & ! strain increment considered significant (used by crystallite to determine whether strain inc is considered significant) defgradTolerance, & ! deviation of deformation gradient that is still allowed (used by CPFEM to determine outdated ffn1) pert_Fg, & ! strain perturbation for FEM Jacobi @@ -106,7 +108,9 @@ subroutine numerics_init() rTol_crystalliteTemperature = 1.0e-6_pReal rTol_crystalliteStress = 1.0e-6_pReal aTol_crystalliteStress = 1.0e-8_pReal ! residuum is in Lp (hence strain on the order of 1e-8 here) - + integrator = 1 + integratorStiffness = 1 + !* RGC parameters: added <<>> with moderate setting absTol_RGC = 1.0e+4 relTol_RGC = 1.0e-3 @@ -181,6 +185,10 @@ subroutine numerics_init() rTol_crystalliteStress = IO_floatValue(line,positions,2) case ('atol_crystallitestress') aTol_crystalliteStress = IO_floatValue(line,positions,2) + case ('integrator') + integrator = IO_intValue(line,positions,2) + case ('integratorstiffness') + integratorStiffness = IO_intValue(line,positions,2) !* RGC parameters: added <<>> case ('atol_rgc') @@ -242,6 +250,8 @@ subroutine numerics_init() write(6,'(a24,x,e8.1)') 'rTol_crystalliteTemp: ',rTol_crystalliteTemperature write(6,'(a24,x,e8.1)') 'rTol_crystalliteStress: ',rTol_crystalliteStress write(6,'(a24,x,e8.1)') 'aTol_crystalliteStress: ',aTol_crystalliteStress + write(6,'(a24,x,i8)') 'integrator: ',integrator + write(6,'(a24,x,i8)') 'integratorStiffness: ',integratorStiffness write(6,*) write(6,'(a24,x,i8)') 'nHomog: ',nHomog @@ -252,23 +262,23 @@ subroutine numerics_init() write(6,*) !* RGC parameters: added <<>> - write(6,'(a24,x,e8.1)') 'aTol_RGC: ',absTol_RGC - write(6,'(a24,x,e8.1)') 'rTol_RGC: ',relTol_RGC - write(6,'(a24,x,e8.1)') 'aMax_RGC: ',absMax_RGC - write(6,'(a24,x,e8.1)') 'rMax_RGC: ',relMax_RGC - write(6,'(a24,x,e8.1)') 'perturbPenalty_RGC: ',pPert_RGC - write(6,'(a24,x,e8.1)') 'relevantMismatch_RGC: ',xSmoo_RGC - write(6,'(a24,x,e8.1)') 'viscosityrate_RGC: ',viscPower_RGC - write(6,'(a24,x,e8.1)') 'viscositymodulus_RGC: ',viscModus_RGC - write(6,'(a24,x,e8.1)') 'maxrelaxation_RGC: ',maxdRelax_RGC - write(6,'(a24,x,e8.1)') 'maxVolDiscrepancy_RGC:',maxVolDiscr_RGC - write(6,'(a24,x,e8.1)') 'volDiscrepancyMod_RGC:',volDiscrMod_RGC - write(6,'(a24,x,e8.1)') 'discrepancyPower_RGC: ',volDiscrPow_RGC + write(6,'(a24,x,e8.1)') 'aTol_RGC: ',absTol_RGC + write(6,'(a24,x,e8.1)') 'rTol_RGC: ',relTol_RGC + write(6,'(a24,x,e8.1)') 'aMax_RGC: ',absMax_RGC + write(6,'(a24,x,e8.1)') 'rMax_RGC: ',relMax_RGC + write(6,'(a24,x,e8.1)') 'perturbPenalty_RGC: ',pPert_RGC + write(6,'(a24,x,e8.1)') 'relevantMismatch_RGC: ',xSmoo_RGC + write(6,'(a24,x,e8.1)') 'viscosityrate_RGC: ',viscPower_RGC + write(6,'(a24,x,e8.1)') 'viscositymodulus_RGC: ',viscModus_RGC + write(6,'(a24,x,e8.1)') 'maxrelaxation_RGC: ',maxdRelax_RGC + write(6,'(a24,x,e8.1)') 'maxVolDiscrepancy_RGC: ',maxVolDiscr_RGC + write(6,'(a24,x,e8.1)') 'volDiscrepancyMod_RGC: ',volDiscrMod_RGC + write(6,'(a24,x,e8.1)') 'discrepancyPower_RGC: ',volDiscrPow_RGC write(6,*) !* Random seeding parameters: added <<>> - write(6,'(a24,x,i8)') 'fixed_seed: ',fixedSeed + write(6,'(a24,x,i8)') 'fixed_seed: ',fixedSeed write(6,*) ! sanity check @@ -294,6 +304,10 @@ subroutine numerics_init() if (rTol_crystalliteTemperature <= 0.0_pReal) call IO_error(276) !! oops !! if (rTol_crystalliteStress <= 0.0_pReal) call IO_error(270) if (aTol_crystalliteStress <= 0.0_pReal) call IO_error(271) + if (integrator <= 0_pInt .or. integrator >= 6_pInt) & + call IO_error(298) + if (integratorStiffness <= 0_pInt .or. integratorStiffness >= 6_pInt) & + call IO_error(298) !* RGC parameters: added <<>> if (absTol_RGC <= 0.0_pReal) call IO_error(272)