diff --git a/code/IO.f90 b/code/IO.f90 index 5b98111df..52e978115 100644 --- a/code/IO.f90 +++ b/code/IO.f90 @@ -1118,6 +1118,8 @@ endfunction msg = '+ Crystallite responds elastically +' case (650) msg = '+ Polar decomposition failed +' + case (700) + msg = '+ unknown crystal symmetry +' case default msg = '+ Unknown warning number... +' end select diff --git a/code/constitutive.f90 b/code/constitutive.f90 index d1eb01ea6..f4f6b679c 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -18,6 +18,8 @@ MODULE constitutive constitutive_subState0, & ! pointer array to microstructure at start of crystallite inc constitutive_state, & ! pointer array to current microstructure (end of converged time step) constitutive_dotState, & ! pointer array to evolution of current microstructure + constitutive_previousDotState, &! pointer array to previous evolution of current microstructure + constitutive_previousDotState2, &! pointer array to previous evolution of current microstructure constitutive_relevantState ! relevant state values integer(pInt), dimension(:,:,:), allocatable :: constitutive_sizeDotState, & ! size of dotState array constitutive_sizeState, & ! size of state array per grain @@ -112,6 +114,8 @@ subroutine constitutive_init() allocate(constitutive_subState0(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) allocate(constitutive_state(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) allocate(constitutive_dotState(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 @@ -132,6 +136,8 @@ subroutine constitutive_init() allocate(constitutive_state(g,i,e)%p(constitutive_j2_sizeState(myInstance))) allocate(constitutive_relevantState(g,i,e)%p(constitutive_j2_sizeState(myInstance))) allocate(constitutive_dotState(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))) 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) @@ -145,6 +151,8 @@ subroutine constitutive_init() allocate(constitutive_state(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) allocate(constitutive_relevantState(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) allocate(constitutive_dotState(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))) 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) @@ -158,6 +166,8 @@ subroutine constitutive_init() allocate(constitutive_state(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance))) allocate(constitutive_relevantState(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance))) allocate(constitutive_dotState(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))) 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) @@ -171,6 +181,8 @@ subroutine constitutive_init() allocate(constitutive_state(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance))) allocate(constitutive_relevantState(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance))) allocate(constitutive_dotState(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))) 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) @@ -190,6 +202,7 @@ subroutine constitutive_init() constitutive_maxSizeDotState = maxval(constitutive_sizeDotState) constitutive_maxSizePostResults = maxval(constitutive_sizePostResults) + !$OMP CRITICAL (write2out) write(6,*) write(6,*) '<<<+- constitutive init -+>>>' write(6,*) '$Id$' @@ -200,6 +213,7 @@ 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) @@ -207,7 +221,8 @@ subroutine constitutive_init() write(6,'(a32,x,7(i5,x))') 'maxSizeState: ', constitutive_maxSizeState write(6,'(a32,x,7(i5,x))') 'maxSizeDotState: ', constitutive_maxSizeDotState write(6,'(a32,x,7(i5,x))') 'maxSizePostResults: ', constitutive_maxSizePostResults - + call flush(6) + !$OMP END CRITICAL (write2out) return endsubroutine @@ -279,8 +294,7 @@ subroutine constitutive_microstructure(Temperature,Fe,Fp,ipc,ip,el) !* Definition of variables integer(pInt), intent(in) :: ipc,ip,el real(pReal), intent(in) :: Temperature -real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: Fe -real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: Fp +real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: Fe, Fp select case (phase_constitution(material_phase(ipc,ip,el))) @@ -362,39 +376,62 @@ subroutine constitutive_collectDotState(Tstar_v, subTstar0_v, Fe, Fp, Temperatur !* OUTPUT: * !* - constitutive_dotState : evolution of state variable * !********************************************************************* - use prec, only: pReal, pInt - use debug - use mesh, only: mesh_NcpElems, mesh_maxNips - use material, only: phase_constitution, material_phase, homogenization_maxNgrains - use constitutive_j2 - use constitutive_phenopowerlaw - use constitutive_dislotwin - use constitutive_nonlocal - implicit none +use prec, only: pReal, pInt +use debug, only: debug_cumDotStateCalls, & + debug_cumDotStateTicks +use mesh, only: mesh_NcpElems, & + mesh_maxNips +use material, only: phase_constitution, & + material_phase, & + homogenization_maxNgrains +use constitutive_j2 +use constitutive_phenopowerlaw +use constitutive_dislotwin +use constitutive_nonlocal -!* Definition of variables - integer(pInt), intent(in) :: ipc,ip,el - real(pReal), intent(in) :: Temperature, subdt - real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: Fe, Fp - real(pReal), dimension(6), intent(in) :: Tstar_v, subTstar0_v +implicit none - select case (phase_constitution(material_phase(ipc,ip,el))) +!*** input variables +integer(pInt), intent(in) :: ipc, ip, el +real(pReal), intent(in) :: Temperature, & + subdt +real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + Fe, & + Fp +real(pReal), dimension(6), intent(in) :: & + Tstar_v, & + subTstar0_v + +!*** local variables +integer(pLongInt) tick, tock, & + tickrate, & + maxticks + +call system_clock(count=tick,count_rate=tickrate,count_max=maxticks) + +select case (phase_constitution(material_phase(ipc,ip,el))) + + case (constitutive_j2_label) + constitutive_dotState(ipc,ip,el)%p = constitutive_j2_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el) - case (constitutive_j2_label) - constitutive_dotState(ipc,ip,el)%p = constitutive_j2_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el) - - case (constitutive_phenopowerlaw_label) - constitutive_dotState(ipc,ip,el)%p = constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el) - - case (constitutive_dislotwin_label) - constitutive_dotState(ipc,ip,el)%p = constitutive_dislotwin_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el) - - case (constitutive_nonlocal_label) - call constitutive_nonlocal_dotState(constitutive_dotState, Tstar_v, subTstar0_v, Fe, Fp, Temperature, subdt, & - constitutive_state, constitutive_subState0, ipc, ip, el) - - end select - return + case (constitutive_phenopowerlaw_label) + constitutive_dotState(ipc,ip,el)%p = constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el) + + case (constitutive_dislotwin_label) + constitutive_dotState(ipc,ip,el)%p = constitutive_dislotwin_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el) + + case (constitutive_nonlocal_label) + call constitutive_nonlocal_dotState(constitutive_dotState, Tstar_v, subTstar0_v, Fe, Fp, Temperature, subdt, & + constitutive_state, constitutive_subState0, subdt, ipc, ip, el) + +end select + +call system_clock(count=tock,count_rate=tickrate,count_max=maxticks) +debug_cumDotStateCalls = debug_cumDotStateCalls + 1_pInt +debug_cumDotStateTicks = debug_cumDotStateTicks + tock-tick +if (tock < tick) debug_cumDotStateTicks = debug_cumDotStateTicks + maxticks + +return endsubroutine @@ -444,7 +481,7 @@ function constitutive_dotTemperature(Tstar_v,Temperature,ipc,ip,el) endfunction -pure function constitutive_postResults(Tstar_v,subTstar0_v,Temperature,dt,subdt,ipc,ip,el) +function constitutive_postResults(Tstar_v, subTstar0_v, Fe, Fp, Temperature, dt, subdt, ipc, ip, el) !********************************************************************* !* return array of constitutive results * !* INPUT: * @@ -454,8 +491,12 @@ pure function constitutive_postResults(Tstar_v,subTstar0_v,Temperature,dt,subdt, !* - ip : current integration point * !* - el : current element * !********************************************************************* - use prec, only: pReal,pInt - use material, only: phase_constitution,material_phase + use prec, only: pReal,pInt + use mesh, only: mesh_NcpElems, & + mesh_maxNips + use material, only: phase_constitution, & + material_phase, & + homogenization_maxNgrains use constitutive_j2 use constitutive_phenopowerlaw use constitutive_dislotwin @@ -466,6 +507,7 @@ pure function constitutive_postResults(Tstar_v,subTstar0_v,Temperature,dt,subdt, integer(pInt), intent(in) :: ipc,ip,el real(pReal), intent(in) :: dt, Temperature, subdt real(pReal), dimension(6), intent(in) :: Tstar_v, subTstar0_v + real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: Fe, Fp real(pReal), dimension(constitutive_sizePostResults(ipc,ip,el)) :: constitutive_postResults constitutive_postResults = 0.0_pReal @@ -481,8 +523,9 @@ pure function constitutive_postResults(Tstar_v,subTstar0_v,Temperature,dt,subdt, constitutive_postResults = constitutive_dislotwin_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el) case (constitutive_nonlocal_label) - constitutive_postResults = constitutive_nonlocal_postResults(Tstar_v, subTstar0_v, Temperature, dt, subdt, constitutive_state,& - constitutive_subState0, constitutive_dotstate, ipc, ip, el) + constitutive_postResults = constitutive_nonlocal_postResults(Tstar_v, subTstar0_v, Fe, Fp, Temperature, dt, subdt, & + constitutive_state, constitutive_subState0, & + constitutive_dotstate, ipc, ip, el) end select return diff --git a/code/constitutive_nonlocal.f90 b/code/constitutive_nonlocal.f90 index 124d35228..b87475c62 100644 --- a/code/constitutive_nonlocal.f90 +++ b/code/constitutive_nonlocal.f90 @@ -79,8 +79,6 @@ real(pReal), dimension(:,:,:), allocatable :: constitutive_nonlocal_ constitutive_nonlocal_forestProjectionScrew, & ! matrix of forest projections of screw dislocations for each instance constitutive_nonlocal_interactionMatrixSlipSlip ! interaction matrix of the different slip systems for each instance -logical :: periodicBC = .false. - CONTAINS !**************************************** @@ -437,7 +435,11 @@ do i = 1,maxNinstance case( 'rho', & 'delta', & 'rho_edge', & + 'rho_edge_pos', & + 'rho_edge_neg', & 'rho_screw', & + 'rho_screw_pos', & + 'rho_screw_neg', & 'excess_rho', & 'excess_rho_edge', & 'excess_rho_screw', & @@ -457,6 +459,34 @@ do i = 1,maxNinstance 'rho_dot_ann_ath', & 'rho_dot_ann_the', & 'rho_dot_flux', & + 'rho_dot_flux_ep', & + 'rho_dot_flux_en', & + 'rho_dot_flux_sp', & + 'rho_dot_flux_sn', & + 'rho_dot_flux_n1_ep', & + 'rho_dot_flux_n1_en', & + 'rho_dot_flux_n1_sp', & + 'rho_dot_flux_n1_sn', & + 'rho_dot_flux_n2_ep', & + 'rho_dot_flux_n2_en', & + 'rho_dot_flux_n2_sp', & + 'rho_dot_flux_n2_sn', & + 'rho_dot_flux_n3_ep', & + 'rho_dot_flux_n3_en', & + 'rho_dot_flux_n3_sp', & + 'rho_dot_flux_n3_sn', & + 'rho_dot_flux_n4_ep', & + 'rho_dot_flux_n4_en', & + 'rho_dot_flux_n4_sp', & + 'rho_dot_flux_n4_sn', & + 'rho_dot_flux_n5_ep', & + 'rho_dot_flux_n5_en', & + 'rho_dot_flux_n5_sp', & + 'rho_dot_flux_n5_sn', & + 'rho_dot_flux_n6_ep', & + 'rho_dot_flux_n6_en', & + 'rho_dot_flux_n6_sp', & + 'rho_dot_flux_n6_sn', & 'd_upper_edge', & 'd_upper_screw', & 'd_upper_dot_edge', & @@ -782,7 +812,6 @@ real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstan neighboring_rhoScrewExcess,&! screw excess dislocation density of my neighbor neighboring_Nedge, & ! total number of edge excess dislocations in my neighbor neighboring_Nscrew -logical flipConnectingVector myInstance = phase_constitutionInstance(material_phase(g,ip,el)) @@ -810,7 +839,7 @@ forall (s = 1:ns) & rhoForest(s) & = dot_product( (rhoEdgePos + rhoEdgeNeg + rhoEdgeDip), constitutive_nonlocal_forestProjectionEdge(1:ns, s, myInstance) ) & + dot_product( (rhoScrewPos + rhoScrewNeg + rhoScrewDip), constitutive_nonlocal_forestProjectionScrew(1:ns, s, myInstance) ) ! calculation of forest dislocation density as projection of screw and edge dislocations -! if (debugger) write(6,'(a23,3(i3,x),/,12(e10.3,x),/)') 'forest dislocation density at ',g,ip,el, rhoForest +! if (debugger) write(6,'(a30,3(i3,x),/,12(e10.3,x),/)') 'forest dislocation density at ',g,ip,el, rhoForest !*** calculate the threshold shear stress for dislocation slip @@ -835,35 +864,8 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el)) neighboring_el = mesh_ipNeighborhood(1,n,ip,el) neighboring_ip = mesh_ipNeighborhood(2,n,ip,el) - flipConnectingVector = .false. - if ( neighboring_el == 0 .or. neighboring_ip == 0 ) then - if (.not. periodicBC) then - cycle - else - flipConnectingVector = .true. - select case (n) - case (1) - neighboring_el = mesh_ipNeighborhood(1,2,ip,el) - neighboring_ip = mesh_ipNeighborhood(2,2,ip,el) - case (2) - neighboring_el = mesh_ipNeighborhood(1,1,ip,el) - neighboring_ip = mesh_ipNeighborhood(2,1,ip,el) - case (3) - neighboring_el = mesh_ipNeighborhood(1,4,ip,el) - neighboring_ip = mesh_ipNeighborhood(2,4,ip,el) - case (4) - neighboring_el = mesh_ipNeighborhood(1,3,ip,el) - neighboring_ip = mesh_ipNeighborhood(2,3,ip,el) - case (5) - neighboring_el = mesh_ipNeighborhood(1,6,ip,el) - neighboring_ip = mesh_ipNeighborhood(2,6,ip,el) - case (6) - neighboring_el = mesh_ipNeighborhood(1,5,ip,el) - neighboring_ip = mesh_ipNeighborhood(2,5,ip,el) - endselect - endif - endif + if ( neighboring_el == 0 .or. neighboring_ip == 0 ) cycle ! deformation gradients needed for mapping between configurations neighboring_F = math_mul33x33(Fe(:,:,g,neighboring_ip,neighboring_el), Fp(:,:,g,neighboring_ip,neighboring_el)) @@ -872,13 +874,8 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el)) neighboring_invFe = math_inv3x3(Fe(:,:,g,neighboring_ip,neighboring_el)) ! calculate connection vector between me and my neighbor in its lattice configuration - if (flipConnectingVector) then - connectingVector = math_mul33x3(neighboring_invFe, math_mul33x3(Favg, & - -(mesh_ipCenterOfGravity(:,neighboring_ip,neighboring_el) - mesh_ipCenterOfGravity(:,ip,el)) ) ) - else - connectingVector = math_mul33x3(neighboring_invFe, math_mul33x3(Favg, & - (mesh_ipCenterOfGravity(:,neighboring_ip,neighboring_el) - mesh_ipCenterOfGravity(:,ip,el)) ) ) - endif + connectingVector = math_mul33x3(neighboring_invFe, math_mul33x3(Favg, & + (mesh_ipCenterOfGravity(:,neighboring_ip,neighboring_el) - mesh_ipCenterOfGravity(:,ip,el)) ) ) ! neighboring dislocation densities neighboring_rhoEdgePos = state(1, neighboring_ip, neighboring_el)%p( 1: ns) @@ -931,14 +928,6 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el)) * math_mul33x33(neighboringSlip2myLattice, math_mul33x33(sigma, transpose(neighboringSlip2myLattice))) ) enddo enddo -! if (debugger) then - ! !$OMP CRITICAL (write2out) - ! write(6,*) '::: constitutive_nonlocal_microstructure at ',g,ip,el - ! write(6,*) - ! write(6,'(a,/,6(f12.5,x),/)') 'Tdislocation_v / MPa', Tdislocation_v/1e6_pReal - ! !$OMP CRITICAL (write2out) -! endif - !********************************************************************** @@ -1071,7 +1060,7 @@ if (debugger) then ! write(6,'(a,/,12(f12.5,x),/)') 'v', v ! write(6,'(a,/,12(f12.5,x),/)') 'gdot total /1e-3',gdotSlip*1e3_pReal ! write(6,'(a,/,3(3(f12.7,x)/))') 'Lp',Lp - call flush(6) + ! call flush(6) !$OMPEND CRITICAL (write2out) endif @@ -1082,7 +1071,8 @@ endsubroutine !********************************************************************* !* rate of change of microstructure * !********************************************************************* -subroutine constitutive_nonlocal_dotState(dotState, Tstar_v, subTstar0_v, Fe, Fp, Temperature, subdt, state, subState0, g,ip,el) +subroutine constitutive_nonlocal_dotState(dotState, Tstar_v, previousTstar_v, Fe, Fp, Temperature, dt_previous, & + state, previousState, timestep, g,ip,el) use prec, only: pReal, & pInt, & @@ -1095,6 +1085,7 @@ use math, only: math_norm3, & math_mul33x33, & math_inv3x3, & math_det3x3, & + math_Mandel6to33, & pi use mesh, only: mesh_NcpElems, & mesh_maxNips, & @@ -1121,15 +1112,16 @@ integer(pInt), intent(in) :: g, & ! current ip, & ! current integration point el ! current element number real(pReal), intent(in) :: Temperature, & ! temperature - subdt ! substepped crystallite time increment + timestep, & ! substepped crystallite time increment + dt_previous ! time increment between previous and current state real(pReal), dimension(6), intent(in) :: Tstar_v, & ! current 2nd Piola-Kirchhoff stress in Mandel notation - subTstar0_v ! 2nd Piola-Kirchhoff stress in Mandel notation at start of crystallite increment + previousTstar_v ! previous 2nd Piola-Kirchhoff stress in Mandel notation real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & Fe, & ! elastic deformation gradient Fp ! plastic deformation gradient type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & state, & ! current microstructural state - subState0 ! microstructural state at start of crystallite increment + previousState ! previous microstructural state !*** input/output variables type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(inout) :: & @@ -1149,7 +1141,8 @@ integer(pInt) myInstance, & ! current s, & ! index of my current slip system sLattice ! index of my current slip system according to lattice order real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4) :: & - rho, & ! dislocation densities (positive/negative screw and edge without dipoles) + rho, & ! current dislocation densities (positive/negative screw and edge without dipoles) + previousRho, & ! previous dislocation densities (positive/negative screw and edge without dipoles) totalRhoDot, & ! total rate of change of dislocation densities thisRhoDot, & ! rate of change of dislocation densities for this mechanism gdot, & ! shear rates @@ -1158,17 +1151,18 @@ real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstan rhoForest, & ! forest dislocation density tauSlipThreshold, & ! threshold shear stress tauSlip, & ! current resolved shear stress - subTauSlip0, & ! resolved shear stress at start of crystallite increment + previousTauSlip, & ! previous resolved shear stress v, & ! dislocation velocity invLambda, & ! inverse of mean free path for dislocations vClimb ! climb velocity of edge dipoles real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),2) :: & - rhoDip, & ! dipole dislocation densities (screw and edge dipoles) + rhoDip, & ! current dipole dislocation densities (screw and edge dipoles) + previousRhoDip, & ! previous dipole dislocation densities (screw and edge dipoles) totalRhoDipDot, & ! total rate of change of dipole dislocation densities thisRhoDipDot, & ! rate of change of dipole dislocation densities for this mechanism dLower, & ! minimum stable dipole distance for edges and screws dUpper, & ! current maximum stable dipole distance for edges and screws - dUpper0, & ! maximum stable dipole distance for edges and screws at start of crystallite increment + previousDUpper, & ! previous maximum stable dipole distance for edges and screws dUpperDot ! rate of change of the maximum stable dipole distance for edges and screws real(pReal), dimension(3,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4) :: & m ! direction of dislocation motion @@ -1176,25 +1170,25 @@ real(pReal), dimension(3,3) :: F, & ! total de neighboring_F, & ! total deformation gradient of my neighbor Favg ! average total deformation gradient of me and my neighbor real(pReal), dimension(6) :: Tdislocation_v, & ! current dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress - subTdislocation0_v ! dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress at start of crystallite increment -real(pReal), dimension(3) :: surfaceNormal ! surface normal of the current interface -real(pReal) norm_surfaceNormal, & ! euclidic norm of the surface normal - area, & ! area of the current interface + previousTdislocation_v ! previous dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress +real(pReal), dimension(3) :: surfaceNormal, & ! surface normal in lattice configuration + surfaceNormal_currentconf ! surface normal in current configuration +real(pReal) area, & ! area of the current interface detFe, & ! determinant of elastic defornmation gradient D ! self diffusion -logical flipAreaNormal + myInstance = phase_constitutionInstance(material_phase(g,ip,el)) myStructure = constitutive_nonlocal_structure(myInstance) ns = constitutive_nonlocal_totalNslip(myInstance) tauSlip = 0.0_pReal -subTauSlip0 = 0.0_pReal +previousTauSlip = 0.0_pReal v = 0.0_pReal gdot = 0.0_pReal dLower = 0.0_pReal dUpper = 0.0_pReal -dUpper0 = 0.0_pReal +previousDUpper = 0.0_pReal dUpperDot = 0.0_pReal totalRhoDot = 0.0_pReal thisRhoDot = 0.0_pReal @@ -1204,11 +1198,13 @@ thisRhoDipDot = 0.0_pReal !*** shortcut to state variables forall (t = 1:4) rho(:,t) = state(g,ip,el)%p((t-1)*ns+1:t*ns) +forall (t = 1:4) previousRho(:,t) = previousState(g,ip,el)%p((t-1)*ns+1:t*ns) forall (c = 1:2) rhoDip(:,c) = state(g,ip,el)%p((3+c)*ns+1:(4+c)*ns) +forall (c = 1:2) previousRhoDip(:,c) = previousState(g,ip,el)%p((3+c)*ns+1:(4+c)*ns) rhoForest = state(g,ip,el)%p(6*ns+1:7*ns) tauSlipThreshold = state(g,ip,el)%p(7*ns+1:8*ns) Tdislocation_v = state(g,ip,el)%p(8*ns+1:8*ns+6) -subTdislocation0_v = subState0(g,ip,el)%p(8*ns+1:8*ns+6) +previousTdislocation_v = previousState(g,ip,el)%p(8*ns+1:8*ns+6) !**************************************************************************** @@ -1218,7 +1214,7 @@ do s = 1,ns ! loop over slip systems sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance) tauSlip(s) = math_mul6x6( Tstar_v + Tdislocation_v, lattice_Sslip_v(:,sLattice,myStructure) ) - subTauSlip0(s) = math_mul6x6( subTstar0_v + subTdislocation0_v, lattice_Sslip_v(:,sLattice,myStructure) ) + previousTauSlip(s) = math_mul6x6( previousTstar_v + previousTdislocation_v, lattice_Sslip_v(:,sLattice,myStructure) ) enddo @@ -1233,13 +1229,15 @@ if (debugger) then !$OMP CRITICAL (write2out) write(6,*) '::: constitutive_nonlocal_dotState at ',g,ip,el write(6,*) - write(6,'(a,/,12(f12.5,x),/)') 'tauSlip / MPa', tauSlip/1e6_pReal - write(6,'(a,/,12(f12.5,x),/)') 'tauSlipThreshold / MPa', tauSlipThreshold/1e6_pReal - write(6,'(a,/,12(e12.3,x),/)') 'v', v - write(6,'(a,/,4(12(f12.5,x),/))') 'gdot / 1e-3', gdot*1e3_pReal - write(6,'(a,/,(12(f12.5,x),/))') 'gdot total/ 1e-3', sum(gdot,2)*1e3_pReal - call flush(6) - !$OMP CRITICAL (write2out) + write(6,'(a,/,3(3(f12.6,x)/))') 'Tdislocation / MPa', math_Mandel6to33(Tdislocation_v/1e6) + write(6,'(a,/,3(3(f12.6,x)/))') 'Tstar / MPa', math_Mandel6to33(Tstar_v/1e6) + write(6,'(a,/,12(f12.5,x),/)') 'tauSlip / MPa', tauSlip / 1e6_pReal + write(6,'(a,/,12(f12.5,x),/)') 'tauSlipThreshold / MPa', tauSlipThreshold / 1e6_pReal + write(6,'(a,/,12(e12.5,x),/)') 'v', v + write(6,'(a,/,4(12(e12.5,x),/))') 'gdot / 1e-3', gdot*1e3_pReal + ! write(6,'(a,/,(12(f12.5,x),/))') 'gdot total/ 1e-3', sum(gdot,2)*1e3_pReal + ! call flush(6) + !$OMPEND CRITICAL (write2out) endif @@ -1252,28 +1250,31 @@ dUpper(:,2) = min( constitutive_nonlocal_Gmod(myInstance) * constitutive_nonloca / ( 8.0_pReal * pi * abs(tauSlip) ), & 1.0_pReal / sqrt( sum(rho,2)+sum(rhoDip,2) ) ) dUpper(:,1) = dUpper(:,2) / ( 1.0_pReal - constitutive_nonlocal_nu(myInstance) ) -dUpper0(:,2) = min( constitutive_nonlocal_Gmod(myInstance) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) & - / ( 8.0_pReal * pi * abs(subTauSlip0) ), & - 1.0_pReal / sqrt( sum(rho,2)+sum(rhoDip,2) ) ) -dUpper0(:,1) = dUpper0(:,2) / ( 1.0_pReal - constitutive_nonlocal_nu(myInstance) ) +previousDUpper(:,2) = min( constitutive_nonlocal_Gmod(myInstance) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) & + / ( 8.0_pReal * pi * abs(previousTauSlip) ), & + 1.0_pReal / sqrt( sum(previousRho,2) + sum(previousRhoDip,2) ) ) +previousDUpper(:,1) = previousDUpper(:,2) / ( 1.0_pReal - constitutive_nonlocal_nu(myInstance) ) -if (subdt > 0) dUpperDot = (dUpper - dUpper0) / subdt +if (dt_previous > 0.0_pReal) dUpperDot = (dUpper - previousDUpper) / dt_previous if (debugger) then !$OMP CRITICAL (write2out) - write(6,'(a,/,2(12(e12.5,x),/))') 'dUpper:',dUpper - write(6,'(a,/,2(12(e12.5,x),/))') 'dUpperDot:',dUpperDot - call flush(6) - !$OMP CRITICAL (write2out) + ! write(6,'(a,/,2(12(f12.5,x),/))') 'dUpper / micron', dUpper*1e6_pReal + ! write(6,'(a,/,2(12(f12.5,x),/))') 'dUpperDot / micron/s', dUpperDot * 1e6_pReal + ! call flush(6) + !$OMPEND CRITICAL (write2out) endif !**************************************************************************** !*** calculate dislocation multiplication -invLambda = sqrt(rhoForest) / constitutive_nonlocal_lambda0PerSlipSystem(:,myInstance) - -thisRhoDot = spread(0.25_pReal * sum(abs(gdot),2) * invLambda / constitutive_nonlocal_burgersPerSlipSystem(:,myInstance), 2, 4) +thisRhoDot(:,1:2) = spread(0.5_pReal * sum(abs(gdot(:,3:4)),2) * sqrt(rhoForest) & + / constitutive_nonlocal_lambda0PerSlipSystem(:,myInstance) & + / constitutive_nonlocal_burgersPerSlipSystem(:,myInstance), 2, 2) +thisRhoDot(:,3:4) = spread(0.5_pReal * sum(abs(gdot(:,1:2)),2) * sqrt(rhoForest) & + / constitutive_nonlocal_lambda0PerSlipSystem(:,myInstance) & + / constitutive_nonlocal_burgersPerSlipSystem(:,myInstance), 2, 2) thisRhoDipDot = 0.0_pReal ! dipoles don't multiplicate totalRhoDot = totalRhoDot + thisRhoDot @@ -1281,9 +1282,8 @@ totalRhoDipDot = totalRhoDipDot + thisRhoDipDot if (debugger) then !$OMP CRITICAL (write2out) - write(6,'(a,/,6(12(e12.5,x),/))') 'dislocation multiplication', thisRhoDot * subdt, thisRhoDipDot * subdt - call flush(6) - !$OMP CRITICAL (write2out) + write(6,'(a,/,4(12(e12.5,x),/))') 'dislocation multiplication', thisRhoDot * timestep + !$OMPEND CRITICAL (write2out) endif @@ -1305,32 +1305,8 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el)) neighboring_el = mesh_ipNeighborhood(1,n,ip,el) neighboring_ip = mesh_ipNeighborhood(2,n,ip,el) - flipAreaNormal = .false. ! if neighbor exists, total deformation gradient is averaged over me and my neighbor - if (periodicBC .and. (neighboring_el == 0 .or. neighboring_ip == 0) ) then - flipAreaNormal = .true. - select case (n) - case (1) - neighboring_el = mesh_ipNeighborhood(1,2,ip,el) - neighboring_ip = mesh_ipNeighborhood(2,2,ip,el) - case (2) - neighboring_el = mesh_ipNeighborhood(1,1,ip,el) - neighboring_ip = mesh_ipNeighborhood(2,1,ip,el) - case (3) - neighboring_el = mesh_ipNeighborhood(1,4,ip,el) - neighboring_ip = mesh_ipNeighborhood(2,4,ip,el) - case (4) - neighboring_el = mesh_ipNeighborhood(1,3,ip,el) - neighboring_ip = mesh_ipNeighborhood(2,3,ip,el) - case (5) - neighboring_el = mesh_ipNeighborhood(1,6,ip,el) - neighboring_ip = mesh_ipNeighborhood(2,6,ip,el) - case (6) - neighboring_el = mesh_ipNeighborhood(1,5,ip,el) - neighboring_ip = mesh_ipNeighborhood(2,5,ip,el) - endselect - endif if ( neighboring_el > 0 .and. neighboring_ip > 0 ) then neighboring_F = math_mul33x33(Fe(:,:,g,neighboring_ip,neighboring_el), Fp(:,:,g,neighboring_ip,neighboring_el)) Favg = 0.5_pReal * (F + neighboring_F) @@ -1338,35 +1314,34 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el)) Favg = F endif - ! calculate the area and the surface normal (of unit length) of the interface in lattice configuration - if (flipAreaNormal) then - surfaceNormal = math_det3x3(Favg) / detFe & - * math_mul33x3(transpose(Fe(:,:,g,ip,el)), math_mul33x3(Favg,-mesh_ipAreaNormal(:,n,ip,el))) - else - surfaceNormal = math_det3x3(Favg) / detFe & - * math_mul33x3(transpose(Fe(:,:,g,ip,el)), math_mul33x3(Favg,mesh_ipAreaNormal(:,n,ip,el))) - endif - norm_surfaceNormal = math_norm3(surfaceNormal) - surfaceNormal = surfaceNormal / norm_surfaceNormal - area = mesh_ipArea(n,ip,el) * norm_surfaceNormal - - lineLength = 0.0_pReal + ! calculate the normal of the interface in current and lattice configuration + surfaceNormal_currentconf = math_det3x3(Favg) * math_mul33x3(math_inv3x3(transpose(Favg)), mesh_ipAreaNormal(:,n,ip,el)) + surfaceNormal = math_mul33x3(transpose(Fe(:,:,g,ip,el)), surfaceNormal_currentconf) / detFe + ! calculate the area of the interface + area = mesh_ipArea(n,ip,el) * math_norm3(surfaceNormal) + + ! normalize the surface normal to unit length + surfaceNormal = surfaceNormal / math_norm3(surfaceNormal) + + lineLength = 0.0_pReal do s = 1,ns ! loop over slip systems do t = 1,4 ! loop over dislocation types if ( sign(1.0_pReal,math_mul3x3(m(:,s,t),surfaceNormal)) == sign(1.0_pReal,gdot(s,t)) ) then lineLength(s,t) = gdot(s,t) / constitutive_nonlocal_burgersPerSlipSystem(s,myInstance) & - * math_mul3x3(m(:,s,t),surfaceNormal) * area ! dislocation line length that leaves this interface per second + * math_mul3x3(m(:,s,t), surfaceNormal) * area & + * constitutive_nonlocal_transmissivity(Fe, ip, el, neighboring_ip, neighboring_el) ! dislocation line length that leaves this interface per second; weighted by a transmissivity factor thisRhoDot(s,t) = thisRhoDot(s,t) - lineLength(s,t) / mesh_ipVolume(ip,el) ! subtract dislocation density rate (= line length over volume) that leaves through an interface from my dotState ... if ( neighboring_el > 0 .and. neighboring_ip > 0 ) then -!***************************************************************************************************** -!*** OMP locking for this neighbor missing -!***************************************************************************************************** - dotState(1,neighboring_ip,neighboring_el)%p((t-1)*ns+s) = dotState(1,neighboring_ip,neighboring_el)%p((t-1)*ns+s) & + !$OMP CRITICAL (dotStateLock) + !$OMP FLUSH (dotState) + dotState(1,neighboring_ip,neighboring_el)%p((t-1)*ns+s) = dotState(1,neighboring_ip,neighboring_el)%p((t-1)*ns+s) & + lineLength(s,t) / mesh_ipVolume(neighboring_ip,neighboring_el) ! ... and add it to the neighboring dotState (if neighbor exists) + !$OMP FLUSH (dotState) + !$OMPEND CRITICAL (dotStateLock) endif endif enddo @@ -1378,9 +1353,8 @@ totalRhoDipDot = totalRhoDipDot + thisRhoDipDot if (debugger) then !$OMP CRITICAL (write2out) - write(6,'(a,/,6(12(e12.5,x),/))') 'dislocation flux', thisRhoDot * subdt, thisRhoDipDot * subdt - call flush(6) - !$OMP CRITICAL (write2out) + write(6,'(a,/,4(12(e12.5,x),/))') 'dislocation flux', thisRhoDot * timestep + !$OMPEND CRITICAL (write2out) endif @@ -1401,9 +1375,8 @@ totalRhoDipDot = totalRhoDipDot + thisRhoDipDot if (debugger) then !$OMP CRITICAL (write2out) - write(6,'(a,/,6(12(e12.5,x),/))') 'dipole formation by glide', thisRhoDot * subdt, thisRhoDipDot * subdt - call flush(6) - !$OMP CRITICAL (write2out) + write(6,'(a,/,6(12(e12.5,x),/))') 'dipole formation by glide', thisRhoDot * timestep, thisRhoDipDot * timestep + !$OMPEND CRITICAL (write2out) endif @@ -1420,9 +1393,8 @@ totalRhoDipDot = totalRhoDipDot + thisRhoDipDot if (debugger) then !$OMP CRITICAL (write2out) - write(6,'(a,/,6(12(e12.5,x),/))') 'athermal dipole annihilation', thisRhoDot * subdt, thisRhoDipDot * subdt - call flush(6) - !$OMP CRITICAL (write2out) + write(6,'(a,/,2(12(e12.5,x),/))') 'athermal dipole annihilation', thisRhoDipDot * timestep + !$OMPEND CRITICAL (write2out) endif @@ -1443,21 +1415,20 @@ totalRhoDipDot = totalRhoDipDot + thisRhoDipDot if (debugger) then !$OMP CRITICAL (write2out) - write(6,'(a,/,6(12(e12.5,x),/))') 'thermally activated dipole annihilation', thisRhoDot * subdt, thisRhoDipDot * subdt - call flush(6) - !$OMP CRITICAL (write2out) + write(6,'(a,/,2(12(e12.5,x),/))') 'thermally activated dipole annihilation', thisRhoDipDot * timestep + !$OMPEND CRITICAL (write2out) endif -!*** formation by stress change = alteration in dUpper +!*** formation/dissociation by stress change = alteration in dUpper thisRhoDipDot = 0.0_pReal -forall (c=1:2, s=1:ns, dUpperDot(s,c) > 0) & ! stress decrease - thisRhoDipDot(s,c) = 8.0_pReal * rho(s,2*c-1) * rho(s,2*c) * dUpper(s,c) * dUpperDot(s,c) -forall (c=1:2, s=1:ns, dUpperDot(s,c) < 0) & ! increased stress - thisRhoDipDot(s,c) = rhoDip(s,c) * dUpperDot(s,c) / (dUpper(s,c) - dLower(s,c)) +! forall (c=1:2, s=1:ns, dUpperDot(s,c) > 0.0_pReal) & ! stress decrease => dipole formation + ! thisRhoDipDot(s,c) = 8.0_pReal * rho(s,2*c-1) * rho(s,2*c) * previousDUpper(s,c) * dUpperDot(s,c) +forall (c=1:2, s=1:ns, dUpperDot(s,c) < 0.0_pReal) & ! increased stress => dipole dissociation + thisRhoDipDot(s,c) = rhoDip(s,c) * dUpperDot(s,c) / (previousDUpper(s,c) - dLower(s,c)) forall (t=1:4) & thisRhoDot(:,t) = -0.5_pReal * thisRhoDipDot(:,(t-1)/2+1) @@ -1467,30 +1438,121 @@ totalRhoDipDot = totalRhoDipDot + thisRhoDipDot if (debugger) then !$OMP CRITICAL (write2out) - write(6,'(a,/,6(12(e12.5,x),/))') 'dipole stability by stress change', thisRhoDot * subdt, thisRhoDipDot * subdt - call flush(6) - !$OMP CRITICAL (write2out) + write(6,'(a,/,6(12(e12.5,x),/))') 'dipole stability by stress change', thisRhoDot * timestep, thisRhoDipDot * timestep + !$OMPEND CRITICAL (write2out) endif !**************************************************************************** !*** assign the rates of dislocation densities to my dotState -dotState(1,ip,el)%p(1:4*ns) = reshape(totalRhoDot,(/4*ns/)) ! one-dimension only (linear list) -dotState(1,ip,el)%p(4*ns+1:6*ns) = reshape(totalRhoDipDot,(/2*ns/)) +dotState(1,ip,el)%p(1:4*ns) = dotState(1,ip,el)%p(1:4*ns) + reshape(totalRhoDot,(/4*ns/)) +dotState(1,ip,el)%p(4*ns+1:6*ns) = dotState(1,ip,el)%p(4*ns+1:6*ns) + reshape(totalRhoDipDot,(/2*ns/)) if (debugger) then !$OMP CRITICAL (write2out) - write(6,'(a,/,4(12(e12.5,x),/))') 'deltaRho:', totalRhoDot * subdt - write(6,'(a,/,2(12(e12.5,x),/))') 'deltaRhoDip:', totalRhoDipDot * subdt - call flush(6) - !$OMP CRITICAL (write2out) + write(6,'(a,/,4(12(e12.5,x),/))') 'deltaRho:', totalRhoDot * timestep + write(6,'(a,/,2(12(e12.5,x),/))') 'deltaRhoDip:', totalRhoDipDot * timestep + !$OMPEND CRITICAL (write2out) endif endsubroutine +!********************************************************************* +!* transmissivity of IP interface * +!********************************************************************* +function constitutive_nonlocal_transmissivity(Fe, ip,el, neighboring_ip,neighboring_el) + +use prec, only: pReal, & + pInt, & + p_vec +use mesh, only: mesh_NcpElems, & + mesh_maxNips +use material, only: homogenization_maxNgrains, & + material_phase, & + phase_constitutionInstance +use math, only: math_inv3x3, & + math_mul33x33, & + math_pDecomposition, & + math_misorientation +use IO, only: IO_warning +implicit none + +!* input variables +real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + Fe ! elastic deformation gradient +integer(pInt), intent(in) :: ip, & ! my IP + el, & ! my element number + neighboring_ip, & ! neighboring IP + neighboring_el ! neighboring elment number + +!* output variables +real(pReal) constitutive_nonlocal_transmissivity ! factor determining the transmissivity of an IP interface for dislocations + +!* local variables +real(pReal), dimension(3,3) :: U, R, & ! polar decomposition of my Fe + neighboring_U, neighboring_R, & ! polar decomposition of my neighbor's Fe + dR ! net misorientation +real(pReal), dimension(3) :: axis ! rotation axis +real(pReal) angle ! misorientation angle +integer(pInt) myInstance, & ! current instance of this constitution + symmetryType +logical error + +! initialize with perfect transmissivity +constitutive_nonlocal_transmissivity = 1.0_pReal + +! only change transmissivity when we have a valid neighbor +if ((neighboring_el > 0) .and. (neighboring_ip > 0)) then + + ! calculate orientations from elastic deformation gradients +! write(6,'(a,/,3(3(f12.8),/))') 'my Fe', Fe(:,:,1,ip,el) +! write(6,'(a,/,3(3(f12.8),/))') 'neighboring Fe', Fe(:,:,1,neighboring_ip,neighboring_el) + call math_pDecomposition(Fe(:,:,1,ip,el), U, R, error) ! polar decomposition of Fe + if (.not. error) & + call math_pDecomposition(Fe(:,:,1,neighboring_ip,neighboring_el), neighboring_U, neighboring_R, error) ! polar decomposition of my neighbors Fe + if (error) then ! polar decomposition failed? + call IO_warning(650, el, ip, 1) + + else + ! choose symmetry type of my crystal structure + myInstance = phase_constitutionInstance(material_phase(1,ip,el)) + select case (constitutive_nonlocal_structureName(myInstance)) + case ('fcc','bcc') + symmetryType = 1_pInt + case ('hex') + symmetryType = 2_pInt + case default + symmetryType = 0_pInt + end select + + ! calculate misorientation + R = transpose(R) ! transpose defines orientation (rotation from current into lattice conf) + neighboring_R = transpose(neighboring_R) + call math_misorientation(axis, angle, dR, R, neighboring_R, symmetryType) +! write(6,'(a,2(x,i3),/,3(3(f12.8),/))') 'my R at',ip,el, R +! write(6,'(a,2(x,i3),/,3(3(f12.8),/))') 'neighboring R at',neighboring_ip, neighboring_el, neighboring_R +! write(6,'(a,x,i3,x,i3,x,f10.3,x,3(f8.5,x),/)') 'constitutive_nonlocal_transmissivity', ip,el, angle, axis + + ! transmissivity depends on misorientation angle + if (angle < 3.0_pReal) then + constitutive_nonlocal_transmissivity = 1.0_pReal + elseif (angle < 10.0_pReal) then + constitutive_nonlocal_transmissivity = 0.5_pReal + else + constitutive_nonlocal_transmissivity = 0.1_pReal + endif + + endif + +endif + +endfunction + + + !********************************************************************* !* rate of change of temperature * !********************************************************************* @@ -1527,73 +1589,112 @@ endfunction !********************************************************************* !* return array of constitutive results * !********************************************************************* -pure function constitutive_nonlocal_postResults(Tstar_v, subTstar0_v, Temperature, dt, subdt, state, subState0, dotState, g, ip, el) +function constitutive_nonlocal_postResults(Tstar_v, previousTstar_v, Fe, Fp, Temperature, dt, dt_previous, & + state, previousState, dotState, g,ip,el) use prec, only: pReal, & pInt, & p_vec -use math, only: math_mul6x6, & +use math, only: math_norm3, & + math_mul6x6, & + math_mul3x3, & + math_mul33x3, & + math_mul33x33, & + math_inv3x3, & + math_det3x3, & + math_Mandel6to33, & pi use mesh, only: mesh_NcpElems, & - mesh_maxNips + mesh_maxNips, & + mesh_element, & + FE_NipNeighbors, & + mesh_ipNeighborhood, & + mesh_ipVolume, & + mesh_ipArea, & + mesh_ipAreaNormal use material, only: homogenization_maxNgrains, & material_phase, & phase_constitutionInstance, & phase_Noutput -use lattice, only: lattice_Sslip_v, & +use lattice, only: lattice_Sslip, & + lattice_Sslip_v, & + lattice_sd, & + lattice_sn, & + lattice_st, & + lattice_maxNslipFamily, & lattice_NslipSystem implicit none !*** input variables -integer(pInt), intent(in) :: g, & ! current grain number - ip, & ! current integration point - el ! current element number -real(pReal), intent(in) :: dt, & ! time increment - subdt, & ! time increment of crystallite substep - Temperature ! temperature -real(pReal), dimension(6), intent(in) :: Tstar_v, & ! 2nd Piola-Kirchhoff stress in Mandel notation - subTstar0_v ! 2nd Piola-Kirchhoff stress in Mandel notation at start of crystallite inc -type(p_vec), dimension(homogenization_maxNgrains, mesh_maxNips, mesh_NcpElems), intent(in) :: & - state, & ! microstructural state - subState0, & ! microstructural state at start of crystallite inc - dotState ! evolution rate of microstructural state +integer(pInt), intent(in) :: g, & ! current grain number + ip, & ! current integration point + el ! current element number +real(pReal), intent(in) :: Temperature, & ! temperature + dt, & ! time increment + dt_previous ! time increment between previous and current state +real(pReal), dimension(6), intent(in) :: Tstar_v, & ! current 2nd Piola-Kirchhoff stress in Mandel notation + previousTstar_v ! previous 2nd Piola-Kirchhoff stress in Mandel notation +real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + Fe, & ! elastic deformation gradient + Fp ! plastic deformation gradient +type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + state, & ! current microstructural state + previousState, & ! previous microstructural state + dotState ! evolution rate of microstructural state !*** output variables real(pReal), dimension(constitutive_nonlocal_sizePostResults(phase_constitutionInstance(material_phase(g,ip,el)))) :: & constitutive_nonlocal_postResults !*** local variables -integer(pInt) myInstance, & ! current instance of this constitution - myStructure, & ! current lattice structure - ns, & ! short notation for the total number of active slip systems - o, & ! index of current output - s, & ! index of current slip system - sLattice, & ! index of my current slip system according to lattice order - cs, & ! constitutive result index - c, & ! character of dislocation - t ! type of dislocation +integer(pInt) myInstance, & ! current instance of this constitution + myStructure, & ! current lattice structure + ns, & ! short notation for the total number of active slip systems + neighboring_el, & ! element number of my neighbor + neighboring_ip, & ! integration point of my neighbor + c, & ! character of dislocation + cs, & ! constitutive result index + n, & ! index of my current neighbor + o, & ! index of current output + t, & ! type of dislocation + s, & ! index of my current slip system + sLattice ! index of my current slip system according to lattice order +real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),6,4) :: & + fluxes ! outgoing fluxes per slipsystem, neighbor and dislocation type real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4) :: & - rho, & ! dislocation densities (positive/negative screw and edge without dipoles) - rhoDot, & ! evolution rate of dislocation densities (positive/negative screw and edge without dipoles) - gdot ! shear rates + rho, & ! current dislocation densities (positive/negative screw and edge without dipoles) + previousRho, & ! previous dislocation densities (positive/negative screw and edge without dipoles) + rhoDot, & ! evolution rate of dislocation densities (positive/negative screw and edge without dipoles) + gdot, & ! shear rates + lineLength ! dislocation line length leaving the current interface real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: & - rhoForest, & ! forest dislocation density - tauSlipThreshold, & ! threshold shear stress - tauSlip, & ! current resolved shear stress - subTauSlip0, & ! resolved shear stress at start of crystallite increment - v, & ! dislocation velocity - invLambda, & ! inverse of mean free path for dislocations - vClimb ! climb velocity of edge dipoles + rhoForest, & ! forest dislocation density + tauSlipThreshold, & ! threshold shear stress + tauSlip, & ! current resolved shear stress + previousTauSlip, & ! previous resolved shear stress + v, & ! dislocation velocity + invLambda, & ! inverse of mean free path for dislocations + vClimb ! climb velocity of edge dipoles real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),2) :: & - rhoDip, & ! dipole dislocation densities (screw and edge dipoles) - rhoDipDot, & ! evolution rate of dipole dislocation densities (screw and edge dipoles) - dLower, & ! minimum stable dipole distance for edges and screws - dUpper, & ! current maximum stable dipole distance for edges and screws - dUpper0, & ! maximum stable dipole distance for edges and screws at start of crystallite increment - dUpperDot ! rate of change of the maximum stable dipole distance for edges and screws -real(pReal), dimension(6) :: Tdislocation_v, & ! current dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress - subTdislocation0_v ! dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress at start of crystallite increment -real(pReal) D ! self diffusion + rhoDip, & ! current dipole dislocation densities (screw and edge dipoles) + previousRhoDip, & ! previous dipole dislocation densities (screw and edge dipoles) + rhoDipDot, & ! evolution rate of dipole dislocation densities (screw and edge dipoles) + dLower, & ! minimum stable dipole distance for edges and screws + dUpper, & ! current maximum stable dipole distance for edges and screws + previousDUpper, & ! previous maximum stable dipole distance for edges and screws + dUpperDot ! rate of change of the maximum stable dipole distance for edges and screws +real(pReal), dimension(3,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4) :: & + m ! direction of dislocation motion +real(pReal), dimension(3,3) :: F, & ! total deformation gradient + neighboring_F, & ! total deformation gradient of my neighbor + Favg ! average total deformation gradient of me and my neighbor +real(pReal), dimension(6) :: Tdislocation_v, & ! current dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress + previousTdislocation_v ! previous dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress +real(pReal), dimension(3) :: surfaceNormal, & ! surface normal in lattice configuration + surfaceNormal_currentconf ! surface normal in current configuration +real(pReal) area, & ! area of the current interface + detFe, & ! determinant of elastic defornmation gradient + D ! self diffusion myInstance = phase_constitutionInstance(material_phase(g,ip,el)) @@ -1606,11 +1707,13 @@ constitutive_nonlocal_postResults = 0.0_pReal ! short hand notations for state variables forall (t = 1:4) rho(:,t) = state(g,ip,el)%p((t-1)*ns+1:t*ns) +forall (t = 1:4) previousRho(:,t) = previousState(g,ip,el)%p((t-1)*ns+1:t*ns) forall (c = 1:2) rhoDip(:,c) = state(g,ip,el)%p((3+c)*ns+1:(4+c)*ns) +forall (c = 1:2) previousRhoDip(:,c) = previousState(g,ip,el)%p((3+c)*ns+1:(4+c)*ns) rhoForest = state(g,ip,el)%p(6*ns+1:7*ns) tauSlipThreshold = state(g,ip,el)%p(7*ns+1:8*ns) Tdislocation_v = state(g,ip,el)%p(8*ns+1:8*ns+6) -subTdislocation0_v = subState0(g,ip,el)%p(8*ns+1:8*ns+6) +previousTdislocation_v = previousState(g,ip,el)%p(8*ns+1:8*ns+6) forall (t = 1:4) rhoDot(:,t) = dotState(g,ip,el)%p((t-1)*ns+1:t*ns) forall (c = 1:2) rhoDipDot(:,c) = dotState(g,ip,el)%p((3+c)*ns+1:(4+c)*ns) @@ -1619,7 +1722,7 @@ forall (c = 1:2) rhoDipDot(:,c) = dotState(g,ip,el)%p((3+c)*ns+1:(4+c)*ns) do s = 1,ns sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance) tauSlip(s) = math_mul6x6( Tstar_v + Tdislocation_v, lattice_Sslip_v(:,sLattice,myStructure) ) - subTauSlip0(s) = math_mul6x6( subTstar0_v + subTdislocation0_v, lattice_Sslip_v(:,sLattice,myStructure) ) + previousTauSlip(s) = math_mul6x6( previousTstar_v + previousTdislocation_v, lattice_Sslip_v(:,sLattice,myStructure) ) enddo v = constitutive_nonlocal_v0PerSlipSystem(:,myInstance) & @@ -1637,17 +1740,64 @@ dUpper(:,2) = min( constitutive_nonlocal_Gmod(myInstance) * constitutive_nonloca / ( 8.0_pReal * pi * abs(tauSlip) ), & 1.0_pReal / sqrt( sum(rho,2)+sum(rhoDip,2) ) ) dUpper(:,1) = dUpper(:,2) / ( 1.0_pReal - constitutive_nonlocal_nu(myInstance) ) -dUpper0(:,2) = min( constitutive_nonlocal_Gmod(myInstance) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) & - / ( 8.0_pReal * pi * abs(subTauSlip0) ), & - 1.0_pReal / sqrt( sum(rho,2)+sum(rhoDip,2) ) ) -dUpper0(:,1) = dUpper0(:,2) / ( 1.0_pReal - constitutive_nonlocal_nu(myInstance) ) +previousDUpper(:,2) = min( constitutive_nonlocal_Gmod(myInstance) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) & + / ( 8.0_pReal * pi * abs(previousTauSlip) ), & + 1.0_pReal / sqrt( sum(previousRho,2) + sum(previousRhoDip,2) ) ) +previousDUpper(:,1) = previousDUpper(:,2) / ( 1.0_pReal - constitutive_nonlocal_nu(myInstance) ) -if (subdt > 0) then - dUpperDot = (dUpper - dUpper0) / subdt -else +if (dt_previous > 0.0_pReal) then + dUpperDot = (dUpper - previousDUpper) / dt_previous +else dUpperDot = 0.0_pReal endif + +! calculate fluxes +m(:,:,1) = lattice_sd(:, constitutive_nonlocal_slipSystemLattice(:,myInstance), myStructure) +m(:,:,2) = -lattice_sd(:, constitutive_nonlocal_slipSystemLattice(:,myInstance), myStructure) +m(:,:,3) = lattice_st(:, constitutive_nonlocal_slipSystemLattice(:,myInstance), myStructure) +m(:,:,4) = -lattice_st(:, constitutive_nonlocal_slipSystemLattice(:,myInstance), myStructure) + +F = math_mul33x33(Fe(:,:,g,ip,el), Fp(:,:,g,ip,el)) +detFe = math_det3x3(Fe(:,:,g,ip,el)) + +fluxes = 0.0_pReal +do n = 1,FE_NipNeighbors(mesh_element(2,el)) ! loop through my neighbors + + neighboring_el = mesh_ipNeighborhood(1,n,ip,el) + neighboring_ip = mesh_ipNeighborhood(2,n,ip,el) + + ! if neighbor exists, total deformation gradient is averaged over me and my neighbor + if ( neighboring_el > 0 .and. neighboring_ip > 0 ) then + neighboring_F = math_mul33x33(Fe(:,:,g,neighboring_ip,neighboring_el), Fp(:,:,g,neighboring_ip,neighboring_el)) + Favg = 0.5_pReal * (F + neighboring_F) + else + Favg = F + endif + + ! calculate the normal of the interface in current and lattice configuration + surfaceNormal_currentconf = math_det3x3(Favg) * math_mul33x3(math_inv3x3(transpose(Favg)), mesh_ipAreaNormal(:,n,ip,el)) + surfaceNormal = math_mul33x3(transpose(Fe(:,:,g,ip,el)), surfaceNormal_currentconf) / detFe + ! normalize the surface normal to unit length + surfaceNormal = surfaceNormal / math_norm3(surfaceNormal) + + ! calculate the area of the interface + area = mesh_ipArea(n,ip,el) * math_norm3(surfaceNormal) + + do s = 1,ns ! loop over slip systems + do t = 1,4 ! loop over dislocation types + if ( sign(1.0_pReal,math_mul3x3(m(:,s,t),surfaceNormal)) == sign(1.0_pReal,gdot(s,t)) ) & + fluxes(s,n,t) = gdot(s,t) / constitutive_nonlocal_burgersPerSlipSystem(s,myInstance) & + * math_mul3x3(m(:,s,t), surfaceNormal) * area & + * constitutive_nonlocal_transmissivity(Fe, ip, el, neighboring_ip, neighboring_el) & + / mesh_ipVolume(ip,el) + enddo + enddo +enddo + + + + do o = 1,phase_Noutput(material_phase(g,ip,el)) select case(constitutive_nonlocal_output(o,myInstance)) @@ -1663,10 +1813,26 @@ do o = 1,phase_Noutput(material_phase(g,ip,el)) constitutive_nonlocal_postResults(cs+1:cs+ns) = rho(:,1) + rho(:,2) cs = cs + ns + case ('rho_edge_pos') + constitutive_nonlocal_postResults(cs+1:cs+ns) = rho(:,1) + cs = cs + ns + + case ('rho_edge_neg') + constitutive_nonlocal_postResults(cs+1:cs+ns) = rho(:,2) + cs = cs + ns + case ('rho_screw') constitutive_nonlocal_postResults(cs+1:cs+ns) = rho(:,3) + rho(:,4) cs = cs + ns + case ('rho_screw_pos') + constitutive_nonlocal_postResults(cs+1:cs+ns) = rho(:,3) + cs = cs + ns + + case ('rho_screw_neg') + constitutive_nonlocal_postResults(cs+1:cs+ns) = rho(:,4) + cs = cs + ns + case ('excess_rho') constitutive_nonlocal_postResults(cs+1:cs+ns) = rho(:,1) - rho(:,2) + rho(:,3) - rho(:,4) cs = cs + ns @@ -1734,18 +1900,18 @@ do o = 1,phase_Noutput(material_phase(g,ip,el)) 4.0_pReal * dUpper(:,c) / constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) & * ( rho(:,2*c-1)*abs(gdot(:,2*c)) + rho(:,2*c)*abs(gdot(:,2*c-1)) ) enddo - do c=1,2 - forall (s=1:ns, dUpperDot(s,c) > 0) & ! dipole formation by stress decrease - constitutive_nonlocal_postResults(cs+s) = constitutive_nonlocal_postResults(cs+s) + & - 8.0_pReal * rho(s,2*c-1) * rho(s,2*c) * dUpper(s,c) * dUpperDot(s,c) - enddo + ! do c=1,2 + ! forall (s=1:ns, dUpperDot(s,c) > 0.0_pReal) & ! dipole formation by stress decrease + ! constitutive_nonlocal_postResults(cs+s) = constitutive_nonlocal_postResults(cs+s) + & + ! 8.0_pReal * rho(s,2*c-1) * rho(s,2*c) * previousDUpper(s,c) * dUpperDot(s,c) + ! enddo cs = cs + ns case ('rho_dot_dip2sgl') do c=1,2 - forall (s=1:ns, dUpperDot(s,c) < 0) & + forall (s=1:ns, dUpperDot(s,c) < 0.0_pReal) & constitutive_nonlocal_postResults(cs+s) = constitutive_nonlocal_postResults(cs+s) - & - rhoDip(s,c) * dUpperDot(s,c) / (dUpper(s,c) - dLower(s,c)) + rhoDip(s,c) * dUpperDot(s,c) / (previousDUpper(s,c) - dLower(s,c)) enddo cs = cs + ns @@ -1769,9 +1935,121 @@ do o = 1,phase_Noutput(material_phase(g,ip,el)) ! !!! cross-slip of screws missing !!! cs = cs + ns + case ('rho_dot_flux_n1_ep') + constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,1,1) + cs = cs + ns + + case ('rho_dot_flux_n1_en') + constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,1,2) + cs = cs + ns + + case ('rho_dot_flux_n1_sp') + constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,1,3) + cs = cs + ns + + case ('rho_dot_flux_n1_sn') + constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,1,4) + cs = cs + ns + + case ('rho_dot_flux_n2_ep') + constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,2,1) + cs = cs + ns + + case ('rho_dot_flux_n2_en') + constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,2,2) + cs = cs + ns + + case ('rho_dot_flux_n2_sp') + constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,2,3) + cs = cs + ns + + case ('rho_dot_flux_n2_sn') + constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,2,4) + cs = cs + ns + + case ('rho_dot_flux_n3_ep') + constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,3,1) + cs = cs + ns + + case ('rho_dot_flux_n3_en') + constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,3,2) + cs = cs + ns + + case ('rho_dot_flux_n3_sp') + constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,3,3) + cs = cs + ns + + case ('rho_dot_flux_n3_sn') + constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,3,4) + cs = cs + ns + + case ('rho_dot_flux_n4_ep') + constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,4,1) + cs = cs + ns + + case ('rho_dot_flux_n4_en') + constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,4,2) + cs = cs + ns + + case ('rho_dot_flux_n4_sp') + constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,4,3) + cs = cs + ns + + case ('rho_dot_flux_n4_sn') + constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,4,4) + cs = cs + ns + + case ('rho_dot_flux_n5_ep') + constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,5,1) + cs = cs + ns + + case ('rho_dot_flux_n5_en') + constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,5,2) + cs = cs + ns + + case ('rho_dot_flux_n5_sp') + constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,5,3) + cs = cs + ns + + case ('rho_dot_flux_n5_sn') + constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,5,4) + cs = cs + ns + + case ('rho_dot_flux_n6_ep') + constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,6,1) + cs = cs + ns + + case ('rho_dot_flux_n6_en') + constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,6,2) + cs = cs + ns + + case ('rho_dot_flux_n6_sp') + constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,6,3) + cs = cs + ns + + case ('rho_dot_flux_n6_sn') + constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,6,4) + cs = cs + ns + + case ('rho_dot_flux_ep') + constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(fluxes(:,:,1), 2) + cs = cs + ns + + case ('rho_dot_flux_en') + constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(fluxes(:,:,2), 2) + cs = cs + ns + + case ('rho_dot_flux_sp') + constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(fluxes(:,:,3), 2) + cs = cs + ns + + case ('rho_dot_flux_sn') + constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(fluxes(:,:,4), 2) + cs = cs + ns + case ('rho_dot_flux') - ! !!! still has to be implemented !!! - constitutive_nonlocal_postResults(cs+1:cs+ns) = 0.0_pReal + constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(fluxes(:,:,1), 2) + sum(fluxes(:,:,2), 2) & + + sum(fluxes(:,:,3), 2) + sum(fluxes(:,:,4), 2) cs = cs + ns case ('d_upper_edge') diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 5cc2869bd..366087465 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -19,7 +19,7 @@ implicit none ! **************************************************************** ! *** General variables for the crystallite calculation *** ! **************************************************************** -integer(pInt), parameter :: crystallite_Nresults = 14_pInt ! phaseID, volume, Euler angles, def gradient +integer(pInt), parameter :: crystallite_Nresults = 5_pInt ! phaseID, volume, Euler angles, def gradient real(pReal), dimension (:,:,:), allocatable :: crystallite_dt, & ! requested time increment of each grain crystallite_subdt, & ! substepped time increment of each grain @@ -50,6 +50,7 @@ real(pReal), dimension (:,:,:,:,:), allocatable :: crystallite_Fe, & crystallite_P ! 1st Piola-Kirchhoff stress per grain real(pReal), dimension (:,:,:,:,:,:,:), allocatable :: crystallite_dPdF, & ! individual dPdF per grain crystallite_fallbackdPdF ! dPdF fallback for non-converged grains (elastic prediction) +real(pReal) crystallite_statedamper ! damping for state update logical, dimension (:,:,:), allocatable :: crystallite_localConstitution, & ! indicates this grain to have purely local constitutive law crystallite_requested, & ! flag to request crystallite calculation @@ -57,7 +58,7 @@ logical, dimension (:,:,:), allocatable :: crystallite_localConstit crystallite_converged, & ! convergence flag crystallite_stateConverged, & ! flag indicating convergence of state crystallite_temperatureConverged, & ! flag indicating convergence of temperature - crystallite_todo ! requested and ontrack but not converged + crystallite_todo ! requested and ontrack but not converged CONTAINS @@ -247,10 +248,15 @@ subroutine crystallite_stressAndItsTangent(updateJaco) math_Plain3333to99 use FEsolving, only: FEsolving_execElem, & FEsolving_execIP, & - theInc - use mesh, only: mesh_element - use material, only: homogenization_Ngrains + theInc, & + cycleCounter + use mesh, only: mesh_element, & + mesh_NcpElems, & + mesh_maxNips + use material, only: homogenization_Ngrains, & + homogenization_maxNgrains use constitutive, only: constitutive_maxSizeState, & + constitutive_maxSizeDotState, & constitutive_sizeState, & constitutive_sizeDotState, & constitutive_state, & @@ -258,9 +264,12 @@ subroutine crystallite_stressAndItsTangent(updateJaco) constitutive_partionedState0, & constitutive_homogenizedC, & constitutive_dotState, & + constitutive_previousDotState, & + constitutive_previousDotState2, & constitutive_collectDotState, & - constitutive_dotTemperature - + constitutive_dotTemperature, & + constitutive_microstructure + implicit none !*** input variables ***! @@ -272,33 +281,47 @@ subroutine crystallite_stressAndItsTangent(updateJaco) real(pReal) myTemperature ! local copy of the temperature 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 - myF, & ! local copy of the deformation gradient - myFp, & ! local copy of the plastic deformation gradient - myInvFp, & ! local copy of the invert of plastic deformation gradient - myFe, & ! local copy of the elastic deformation gradient - myLp, & ! local copy of the plastic velocity gradient - myP ! local copy of the 1st Piola-Kirchhoff stress tensor - real(pReal), dimension(6) :: myTstar_v ! local copy of the 2nd Piola-Kirchhoff stress vector - real(pReal), dimension(constitutive_maxSizeState) :: myState, & ! local copy of the state - myDotState ! local copy of dotState + 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, & ! element index - i, & ! integration point index - g, & ! grain index + integer(pInt) e, ee, & ! element index + i, ii, & ! integration point index + g, gg, & ! grain index k, & l, & + comp, & myNgrains, & mySizeState, & mySizeDotState + integer(pInt), dimension(2,9) :: 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) :: dPdF_pos,dPdF_neg - + real(pReal), dimension(constitutive_maxSizeDotState) :: delta_dotState1, & ! difference between current and previous dotstate + delta_dotState2 ! difference between previousDotState and previousDotState2 + real(pReal) dot_prod12, & + dot_prod22 + real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & + storedF, & + storedFp, & + storedInvFp, & + storedFe, & + storedLp, & + storedP + real(pReal), dimension(6,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & + storedTstar_v + real(pReal), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & + storedTemperature + real(pReal), dimension(constitutive_maxSizeState,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & + storedState, & + storedDotState + logical, dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & + storedConvergenceFlag + logical, dimension(3,3) :: mask + ! ------ initialize to starting condition ------ !$OMP CRITICAL (write2out) @@ -311,7 +334,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) ! 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 myNgrains = homogenization_Ngrains(mesh_element(3,e)) @@ -341,13 +365,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) NiterationCrystallite = 0_pInt do while (any(crystallite_subStep(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) > subStepMinCryst)) ! cutback loop for crystallites - - if (any(.not. crystallite_converged(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) .and. & ! any non-converged grain - .not. crystallite_localConstitution(:,:,FEsolving_execELem(1):FEsolving_execElem(2))) ) & ! has non-local constitution? - crystallite_converged(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) = & - crystallite_converged(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) .and. & - crystallite_localConstitution(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) ! reset non-local grains' convergence status - + !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed myNgrains = homogenization_Ngrains(mesh_element(3,e)) @@ -357,7 +375,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) if (crystallite_converged(g,i,e)) then if (debugger) then !$OMP CRITICAL (write2out) - write(6,'(a21,f10.8,a33,f10.8,a35)') 'winding forward from ', & + write(6,'(a21,f10.8,a32,f10.8,a35)') 'winding forward from ', & crystallite_subFrac(g,i,e),' to current crystallite_subfrac ', & crystallite_subFrac(g,i,e)+crystallite_subStep(g,i,e),' in crystallite_stressAndItsTangent' write(6,*) @@ -396,8 +414,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) endif endif - crystallite_onTrack(g,i,e) = crystallite_subStep(g,i,e) > subStepMinCryst ! still on track or already done (beyond repair) - if (crystallite_onTrack(g,i,e)) then ! specify task (according to substep) + crystallite_onTrack(g,i,e) = crystallite_subStep(g,i,e) > subStepMinCryst ! still on track or already done (beyond repair) + if (crystallite_onTrack(g,i,e)) then ! specify task (according to substep) 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)) @@ -411,9 +429,11 @@ subroutine crystallite_stressAndItsTangent(updateJaco) !$OMPEND PARALLEL DO crystallite_todo = ( crystallite_requested & - .and. crystallite_onTrack & - .and. .not. crystallite_converged) - + .and. crystallite_onTrack & + .and. .not. crystallite_converged) + + crystallite_statedamper = 1.0_pReal + ! --+>> preguess for state <<+-- ! ! incrementing by crystallite_subdt @@ -427,9 +447,13 @@ 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 - if (crystallite_todo(g,i,e)) & ! all undone crystallites - constitutive_dotState(g,i,e)%p = 0.0_pReal ! zero out dotState - enddo; enddo; enddo + if (crystallite_todo(g,i,e)) then ! all undone crystallites + call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe, crystallite_Fp, 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 @@ -437,30 +461,27 @@ subroutine crystallite_stressAndItsTangent(updateJaco) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains ! debugger = (e == 1 .and. i == 1 .and. g == 1) - if (crystallite_todo(g,i,e)) & ! all undone crystallites + 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_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 ! debugger = (e == 1 .and. i == 1 .and. g == 1) - if (crystallite_todo(g,i,e)) then ! all undone crystallites + 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) = .false. ! force at least one iteration step even if state already converged endif enddo; enddo; enddo !$OMPEND PARALLEL DO - if (debugger) then - !$OMP CRITICAL (write2out) - write(6,*) count(crystallite_onTrack(1,:,:)),'IPs onTrack after preguess for state' - !$OMPEND CRITICAL (write2out) - endif ! --+>> state loop <<+-- @@ -484,7 +505,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains ! debugger = (e == 1 .and. i == 1 .and. g == 1) - if (crystallite_todo(g,i,e)) & ! all undone crystallites + if (crystallite_todo(g,i,e)) & ! all undone crystallites crystallite_onTrack(g,i,e) = crystallite_integrateStress(g,i,e) enddo enddo @@ -492,14 +513,20 @@ subroutine crystallite_stressAndItsTangent(updateJaco) !$OMPEND PARALLEL DO if (debugger) then !$OMP CRITICAL (write2out) - write(6,*) count(crystallite_onTrack(1,:,:)),'IPs onTrack after stress integration' + write(6,*) count(crystallite_onTrack(:,:,:)),'grains onTrack after stress integration' !$OMPEND CRITICAL (write2out) endif - crystallite_todo = crystallite_todo .and. crystallite_onTrack - if (any(.not. crystallite_onTrack .and. .not. crystallite_localConstitution)) & + crystallite_todo = crystallite_todo .and. crystallite_onTrack ! continue with non-broken grains + if (any(.not. crystallite_onTrack .and. .not. crystallite_localConstitution)) & ! any non-local is broken? crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! all nonlocal crystallites can be skipped + if (debugger) then + !$OMP CRITICAL (write2out) + write(6,*) count(crystallite_todo(:,:,:)),'grains todo after stress integration' + !$OMPEND CRITICAL (write2out) + endif + ! --+>> state integration <<+-- ! ! incrementing by crystallite_subdt @@ -513,20 +540,34 @@ 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 - if (crystallite_todo(g,i,e)) & ! all undone crystallites + 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 ! debugger = (e == 1 .and. i == 1 .and. g == 1) - if (crystallite_todo(g,i,e)) & ! all undone crystallites + 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_subdt(g,i,e), g, i, e) + crystallite_subdt(g,i,e), g, i, e) + delta_dotState1 = constitutive_dotState(g,i,e)%p - constitutive_previousDotState(g,i,e)%p + delta_dotState2 = constitutive_previousDotState(g,i,e)%p - constitutive_previousDotState2(g,i,e)%p + dot_prod12 = dot_product(delta_dotState1, delta_dotState2) + dot_prod22 = dot_product(delta_dotState2, delta_dotState2) + 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 = min(crystallite_statedamper, & + 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 @@ -550,16 +591,21 @@ subroutine crystallite_stressAndItsTangent(updateJaco) enddo enddo !$OMPEND PARALLEL DO + + if (debugger) then + !$OMP CRITICAL (write2out) + write(6,*) count(crystallite_converged(:,:,:)),'grains converged after state integration no.', NiterationState + !$OMPEND CRITICAL (write2out) + endif + 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. crystallite_onTrack .and. .not. crystallite_converged - if (any(.not. crystallite_onTrack .and. .not. crystallite_localConstitution)) & - crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! all nonlocal crystallites can be skipped + crystallite_todo = crystallite_todo .and. .not. crystallite_converged ! skip all converged if (debugger) then !$OMP CRITICAL (write2out) - write(6,*) count(crystallite_onTrack(1,:,:)),'IPs onTrack after state update' - write(6,*) count(crystallite_converged(1,:,:)),'IPs converged' - write(6,*) count(crystallite_todo(1,:,:)),'IPs todo' + 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 @@ -570,9 +616,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) enddo ! cutback loop - ! write (6,'(a,/,32(L,x))') 'crystallite_todo',crystallite_todo - ! write (6,'(a,/,32(L,x))') 'crystallite_converged',crystallite_converged - ! ------ check for non-converged crystallites ------ !$OMP PARALLEL DO @@ -588,7 +631,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) math_mul66x6( 0.5_pReal*constitutive_homogenizedC(g,i,e), & math_Mandel33to6( math_mul33x33(transpose(Fe_guess),Fe_guess) - math_I3 ) & ) & - ) + ) crystallite_P(:,:,g,i,e) = math_mul33x33(Fe_guess,math_mul33x33(Tstar,transpose(invFp))) endif enddo @@ -596,168 +639,296 @@ subroutine crystallite_stressAndItsTangent(updateJaco) enddo !$OMPEND PARALLEL DO + ! --+>> stiffness calculation <<+-- + + if(updateJaco) then ! Jacobian required + + crystallite_statedamper = 1.0_pReal + + 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 + 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 + storedState(1:mySizeState,g,i,e) = constitutive_state(g,i,e)%p ! remember unperturbed, converged state, ... + storedDotState(1:mySizeDotState,g,i,e) = constitutive_dotState(g,i,e)%p ! ... dotStates, ... + storedTemperature(g,i,e) = crystallite_Temperature(g,i,e) ! ... Temperature, ... + storedF(:,:,g,i,e) = crystallite_subF(:,:,g,i,e) ! ... and kinematics + storedFp(:,:,g,i,e) = crystallite_Fp(:,:,g,i,e) + storedInvFp(:,:,g,i,e) = crystallite_invFp(:,:,g,i,e) + storedFe(:,:,g,i,e) = crystallite_Fe(:,:,g,i,e) + storedLp(:,:,g,i,e) = crystallite_Lp(:,:,g,i,e) + storedTstar_v(:,g,i,e) = crystallite_Tstar_v(:,g,i,e) + storedP(:,:,g,i,e) = crystallite_P(:,:,g,i,e) + storedConvergenceFlag(g,i,e) = crystallite_converged(g,i,e) + enddo; enddo; enddo - if(updateJaco) then ! Jacobian required - !$OMP PARALLEL DO - do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed + if (all(crystallite_localConstitution) .or. theInc < 2) then ! all grains have local constitution, so local convergence of perturbed grain is sufficient + + !$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_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) then + !$OMP CRITICAL (write2out) + write (6,*) '#############' + write (6,*) 'central solution of cryst_StressAndTangent' + write (6,*) '#############' + write (6,'(a,/,3(3(f12.4,x)/))') ' P of 1 1 1',storedP(1:3,:,g,i,e)/1e6 + write (6,'(a,/,3(3(f12.8,x)/))') ' Fp of 1 1 1',storedFp(1:3,:,g,i,e) + write (6,'(a,/,3(3(f12.8,x)/))') ' Lp of 1 1 1',storedLp(1:3,:,g,i,e) + !$OMPEND CRITICAL (write2out) + endif + + if (pert_method == 1_pInt .or. pert_method == 3_pInt) then ! <<< when forward or central difference is desired >>> + do k = 1,3 ! perturbation... + do l = 1,3 ! ...components to the positive direction + crystallite_subF(k,l,g,i,e) = crystallite_subF(k,l,g,i,e) + pert_Fg ! perturb single component + if (debugger) then + !$OMP CRITICAL (write2out) + write (6,'(i1,x,i1)') k,l + write (6,'(a,/,3(3(f12.6,x)/))') 'pertF of 1 1 1',crystallite_subF(1:3,:,g,i,e) + !$OMPEND CRITICAL (write2out) + endif + onTrack = .true. + converged = .false. + NiterationState = 0_pInt + do while(.not. converged .and. onTrack .and. NiterationState < nState) ! keep cycling until done (potentially non-converged) + NiterationState = NiterationState + 1_pInt + onTrack = crystallite_integrateStress(g,i,e) ! stress of perturbed situation (overwrites _P,_Tstar_v,_Fp,_Lp,_Fe) + 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_subdt(g,i,e), g, i, e) + stateConverged = crystallite_updateState(g,i,e) ! update state + temperatureConverged = crystallite_updateTemperature(g,i,e) ! update temperature + converged = stateConverged .and. temperatureConverged + endif + if (debugger) then + !$OMP CRITICAL (write2out) + write (6,*) '-------------' + write (6,'(a,x,l,x,l)') 'ontrack + converged:',onTrack,converged + write (6,'(a,/,3(3(f12.4,x)/))') 'pertP of 1 1 1',crystallite_P(1:3,:,g,i,e)/1e6 + write (6,'(a,/,3(3(f12.4,x)/))') 'DP of 1 1 1',(crystallite_P(1:3,:,g,i,e)-storedP(1:3,:,g,i,e))/1e6 + !$OMPEND CRITICAL (write2out) + endif + enddo + if (converged) then ! converged state warrants stiffness update + dPdF_pos(:,:,k,l) = (crystallite_P(:,:,g,i,e) - storedP(:,:,g,i,e))/pert_Fg ! tangent dP_ij/dFg_kl + if (pert_method == 1_pInt) crystallite_dPdF(:,:,k,l,g,i,e) = dPdF_pos(:,:,k,l) + endif + do ee = FEsolving_execElem(1),FEsolving_execElem(2) + myNgrains = homogenization_Ngrains(mesh_element(3,ee)) + do ii = FEsolving_execIP(1,ee),FEsolving_execIP(2,ee) + do gg = 1,myNgrains + mySizeState = constitutive_sizeState(gg,ii,ee) ! number of state variables for this grain + mySizeDotState = constitutive_sizeDotState(gg,ii,ee) ! number of dotStates for this grain + constitutive_state(gg,ii,ee)%p = storedState(1:mySizeState,gg,ii,ee) + constitutive_dotState(gg,ii,ee)%p = storedDotState(1:mySizeDotState,gg,ii,ee) + crystallite_Temperature(gg,ii,ee) = storedTemperature(gg,ii,ee) + crystallite_subF(:,:,gg,ii,ee) = storedF(:,:,gg,ii,ee) + crystallite_Fp(:,:,gg,ii,ee) = storedFp(:,:,gg,ii,ee) + crystallite_invFp(:,:,gg,ii,ee) = storedInvFp(:,:,gg,ii,ee) + crystallite_Fe(:,:,gg,ii,ee) = storedFe(:,:,gg,ii,ee) + crystallite_Lp(:,:,gg,ii,ee) = storedLp(:,:,gg,ii,ee) + crystallite_Tstar_v(:,gg,ii,ee) = storedTstar_v(:,gg,ii,ee) + crystallite_P(:,:,gg,ii,ee) = storedP(:,:,gg,ii,ee) + enddo; enddo; enddo + !$OMP CRITICAL (out) + debug_StiffnessStateLoopDistribution(NiterationState) = & + debug_StiffnessstateLoopDistribution(NiterationState) + 1 + !$OMPEND CRITICAL (out) + enddo + enddo + endif + + if (pert_method == 2_pInt .or. pert_method == 3_pInt) then ! <<< when backward or central difference is desired >>> + do k = 1,3 ! perturbation... + do l = 1,3 ! ...components to the negative direction + crystallite_subF(k,l,g,i,e) = crystallite_subF(k,l,g,i,e) - pert_Fg ! perturb single component + if (debugger) then + !$OMP CRITICAL (write2out) + write (6,'(i1,x,i1)') k,l + write (6,'(a,/,3(3(f12.6,x)/))') 'pertF of 1 1 1',crystallite_subF(1:3,:,g,i,e) + !$OMPEND CRITICAL (write2out) + endif + onTrack = .true. + converged = .false. + NiterationState = 0_pInt + do while(.not. converged .and. onTrack .and. NiterationState < nState) ! keep cycling until done (potentially non-converged) + NiterationState = NiterationState + 1_pInt + onTrack = crystallite_integrateStress(g,i,e) ! stress of perturbed situation (overwrites _P,_Tstar_v,_Fp,_Lp,_Fe) + 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_subdt(g,i,e), g, i, e) + stateConverged = crystallite_updateState(g,i,e) ! update state + temperatureConverged = crystallite_updateTemperature(g,i,e) ! update temperature + converged = stateConverged .and. temperatureConverged + endif + if (debugger) then + !$OMP CRITICAL (write2out) + write (6,*) '-------------' + write (6,'(a,x,l,x,l)') 'ontrack + converged:',onTrack,converged + write (6,'(a,/,3(3(f12.4,x)/))') 'pertP of 1 1 1',crystallite_P(1:3,:,g,i,e)/1e6 + write (6,'(a,/,3(3(f12.4,x)/))') 'DP of 1 1 1',(crystallite_P(1:3,:,g,i,e)-storedP(1:3,:,g,i,e))/1e6 + !$OMPEND CRITICAL (write2out) + endif + enddo + if (converged) then ! converged state warrants stiffness update + dPdF_neg(:,:,k,l) = (storedP(:,:,g,i,e) - crystallite_P(:,:,g,i,e))/pert_Fg ! tangent dP_ij/dFg_kl + if (pert_method == 2_pInt) crystallite_dPdF(:,:,k,l,g,i,e) = dPdF_neg(:,:,k,l) + endif + do ee = FEsolving_execElem(1),FEsolving_execElem(2) + myNgrains = homogenization_Ngrains(mesh_element(3,ee)) + do ii = FEsolving_execIP(1,ee),FEsolving_execIP(2,ee) + do gg = 1,myNgrains + mySizeState = constitutive_sizeState(gg,ii,ee) ! number of state variables for this grain + mySizeDotState = constitutive_sizeDotState(gg,ii,ee) ! number of dotStates for this grain + constitutive_state(gg,ii,ee)%p = storedState(1:mySizeState,gg,ii,ee) + constitutive_dotState(gg,ii,ee)%p = storedDotState(1:mySizeDotState,gg,ii,ee) + crystallite_Temperature(gg,ii,ee) = storedTemperature(gg,ii,ee) + crystallite_subF(:,:,gg,ii,ee) = storedF(:,:,gg,ii,ee) + crystallite_Fp(:,:,gg,ii,ee) = storedFp(:,:,gg,ii,ee) + crystallite_invFp(:,:,gg,ii,ee) = storedInvFp(:,:,gg,ii,ee) + crystallite_Fe(:,:,gg,ii,ee) = storedFe(:,:,gg,ii,ee) + crystallite_Lp(:,:,gg,ii,ee) = storedLp(:,:,gg,ii,ee) + crystallite_Tstar_v(:,gg,ii,ee) = storedTstar_v(:,gg,ii,ee) + crystallite_P(:,:,gg,ii,ee) = storedP(:,:,gg,ii,ee) + enddo; enddo; enddo + !$OMP CRITICAL (out) + debug_StiffnessStateLoopDistribution(NiterationState) = & + debug_StiffnessstateLoopDistribution(NiterationState) + 1 + !$OMPEND CRITICAL (out) + enddo + enddo + endif + + if (pert_method == 3_pInt) crystallite_dPdF(:,:,:,:,g,i,e) = 0.5_pReal*(dPdF_neg + dPdF_pos) + + else ! grain did not converge + crystallite_dPdF(:,:,:,:,g,i,e) = crystallite_fallbackdPdF(:,:,:,:,g,i,e) ! use (elastic) fallback + endif ! grain convergence + endif ! grain request + enddo ! grain loop + enddo ! ip loop + enddo ! element loop + !$OMPEND PARALLEL DO + + elseif (any(.not. crystallite_localConstitution)) then ! if any nonlocal grain present, we have to do a full loop over all grains after each perturbance + + 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) ! iterate over IPs of this element to be processed + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do g = 1,myNgrains - 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 - 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 - myState(1:mySizeState) = constitutive_state(g,i,e)%p ! remember unperturbed, converged state, ... - myDotState(1:mySizeDotState) = constitutive_dotState(g,i,e)%p ! ... dotStates, ... - myTemperature = crystallite_Temperature(g,i,e) ! ... Temperature, ... - myF = crystallite_subF(:,:,g,i,e) ! ... and kinematics - myFp = crystallite_Fp(:,:,g,i,e) - myInvFp = crystallite_invFp(:,:,g,i,e) - myFe = crystallite_Fe(:,:,g,i,e) - myLp = crystallite_Lp(:,:,g,i,e) - myTstar_v = crystallite_Tstar_v(:,g,i,e) - myP = crystallite_P(:,:,g,i,e) - if (debugger) then - !$OMP CRITICAL (write2out) - write (6,*) '#############' - write (6,*) 'central solution of cryst_StressAndTangent' - write (6,*) '#############' - write (6,'(a,/,3(3(f12.4,x)/))') ' P of 1 1 1',myP(1:3,:)/1e6 - write (6,'(a,/,3(3(f12.8,x)/))') ' Fp of 1 1 1',myFp(1:3,:) - write (6,'(a,/,3(3(f12.8,x)/))') ' Lp of 1 1 1',myLp(1:3,:) - write (6,'(a,/,16(6(e12.4,x)/),2(f12.4,x))') 'state of 1 1 1',myState/1e6 - !$OMPEND CRITICAL (write2out) - endif + + ! perturb components in the order of biggest change in F (-> component with biggest change in F is perturbed in first cycle, component with second biggest change in next cycle, ...) + mask = .true. + do comp = 1,9 + kl(:,comp) = maxloc(abs(crystallite_subF(:,:,g,i,e)-crystallite_F0(:,:,g,i,e)), mask) + mask(kl(1,comp),kl(2,comp)) = .false. + enddo + k = kl(1,mod((cycleCounter-1)/2+1,9)) + l = kl(2,mod((cycleCounter-1)/2+1,9)) + + crystallite_subF(k,l,g,i,e) = crystallite_subF(k,l,g,i,e) + pert_Fg ! perturb single component + + NiterationState = 0_pInt + crystallite_todo = .true. + do while ( any(crystallite_todo(:,:,FEsolving_execELem(1):FEsolving_execElem(2))) & + .and. NiterationState < nState) + NiterationState = NiterationState + 1_pInt -! begin perturbation of components of F - if (pert_method == 1_pInt .or. pert_method == 3_pInt) then ! <<< when forward or central difference is desired >>> - do k = 1,3 ! perturbation... - do l = 1,3 ! ...components to the positive direction - crystallite_subF(:,:,g,i,e) = myF ! initialize perturbed F to match converged - crystallite_subF(k,l,g,i,e) = crystallite_subF(k,l,g,i,e) + pert_Fg ! perturb single component - if (debugger) then - !$OMP CRITICAL (write2out) - write (6,*) '=============' - write (6,'(i1,x,i1)') k,l - write (6,*) '=============' - write (6,'(a,/,3(3(f12.6,x)/))') 'pertF of 1 1 1',crystallite_subF(1:3,:,g,i,e) - !$OMPEND CRITICAL (write2out) - endif - onTrack = .true. - converged = .false. - NiterationState = 0_pInt - do while(.not. converged .and. onTrack .and. NiterationState < nState) ! keep cycling until done (potentially non-converged) - NiterationState = NiterationState + 1_pInt - onTrack = crystallite_integrateStress(g,i,e) ! stress of perturbed situation (overwrites _P,_Tstar_v,_Fp,_Lp,_Fe) - if (onTrack) 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_subdt(g,i,e), g, i, e) + do ee = FEsolving_execElem(1),FEsolving_execElem(2) + myNgrains = homogenization_Ngrains(mesh_element(3,ee)) + do ii = FEsolving_execIP(1,ee),FEsolving_execIP(2,ee) + do gg = 1,myNgrains + if (crystallite_todo(gg,ii,ee)) & + crystallite_onTrack(gg,ii,ee) = crystallite_integrateStress(gg,ii,ee) ! stress integration + enddo; enddo; enddo + + crystallite_todo = crystallite_todo .and. crystallite_onTrack ! continue with non-broken grains + if (any(.not. crystallite_onTrack .and. .not. crystallite_localConstitution)) & ! any non-local is broken? + crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! all nonlocal crystallites can be skipped + + do ee = FEsolving_execElem(1),FEsolving_execElem(2) + myNgrains = homogenization_Ngrains(mesh_element(3,ee)) + do ii = FEsolving_execIP(1,ee),FEsolving_execIP(2,ee) + do gg = 1,myNgrains + if (crystallite_todo(gg,ii,ee)) & + constitutive_dotState(gg,ii,ee)%p = 0.0_pReal ! zero out dotState + enddo; enddo; enddo - stateConverged = crystallite_updateState(g,i,e) ! update state - temperatureConverged = crystallite_updateTemperature(g,i,e) ! update temperature - converged = stateConverged .and. temperatureConverged - endif - if (debugger) then - !$OMP CRITICAL (write2out) - write (6,*) '-------------' - write (6,'(a,x,l,x,l)') 'ontrack + converged:',onTrack,converged - write (6,'(a,/,3(3(f12.4,x)/))') 'pertP of 1 1 1',crystallite_P(1:3,:,g,i,e)/1e6 - write (6,'(a,/,3(3(f12.4,x)/))') 'DP of 1 1 1',(crystallite_P(1:3,:,g,i,e)-myP(1:3,:))/1e6 - write (6,'(a,/,16(6(e12.4,x)/),/,2(f12.4,x))') 'state of 1 1 1',constitutive_state(g,i,e)%p/1e6 - write (6,'(a,/,16(6(e12.4,x)/),/,2(f12.4,x))') 'Dstate of 1 1 1',(constitutive_state(g,i,e)%p-myState)/1e6 - !$OMPEND CRITICAL (write2out) - endif - enddo - if (converged) then ! converged state warrants stiffness update - dPdF_pos(:,:,k,l) = (crystallite_P(:,:,g,i,e) - myP)/pert_Fg ! tangent dP_ij/dFg_kl - if (pert_method == 1_pInt) crystallite_dPdF(:,:,k,l,g,i,e) = dPdF_pos(:,:,k,l) - endif - constitutive_state(g,i,e)%p = myState ! restore unperturbed, converged state, ... - constitutive_dotState(g,i,e)%p = myDotState ! ... dotState, ... - crystallite_Temperature(g,i,e) = myTemperature ! ... temperature, ... - crystallite_Fp(:,:,g,i,e) = myFp ! ... and kinematics - crystallite_invFp(:,:,g,i,e) = myInvFp - crystallite_Fe(:,:,g,i,e) = myFe - crystallite_Lp(:,:,g,i,e) = myLp - crystallite_Tstar_v(:,g,i,e) = myTstar_v - crystallite_P(:,:,g,i,e) = myP - !$OMP CRITICAL (out) - debug_StiffnessStateLoopDistribution(NiterationState) = & - debug_StiffnessstateLoopDistribution(NiterationState) + 1 - !$OMPEND CRITICAL (out) - enddo - enddo - endif - if (pert_method == 2_pInt .or. pert_method == 3_pInt) then ! <<< when backward or central difference is desired >>> - do k = 1,3 ! perturbation... - do l = 1,3 ! ...components to the negative direction - crystallite_subF(:,:,g,i,e) = myF ! initialize perturbed F to match converged - crystallite_subF(k,l,g,i,e) = crystallite_subF(k,l,g,i,e) - pert_Fg ! perturb single component - if (debugger) then - !$OMP CRITICAL (write2out) - write (6,*) '=============' - write (6,'(i1,x,i1)') k,l - write (6,*) '=============' - write (6,'(a,/,3(3(f12.6,x)/))') 'pertF of 1 1 1',crystallite_subF(1:3,:,g,i,e) - !$OMPEND CRITICAL (write2out) - endif - onTrack = .true. - converged = .false. - NiterationState = 0_pInt - do while(.not. converged .and. onTrack .and. NiterationState < nState) ! keep cycling until done (potentially non-converged) - NiterationState = NiterationState + 1_pInt - onTrack = crystallite_integrateStress(g,i,e) ! stress of perturbed situation (overwrites _P,_Tstar_v,_Fp,_Lp,_Fe) - if (onTrack) 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_subdt(g,i,e), g, i, e) + do ee = FEsolving_execElem(1),FEsolving_execElem(2) + myNgrains = homogenization_Ngrains(mesh_element(3,ee)) + do ii = FEsolving_execIP(1,ee),FEsolving_execIP(2,ee) + do gg = 1,myNgrains + if (crystallite_todo(gg,ii,ee)) & + call constitutive_collectDotState(crystallite_Tstar_v(:,gg,ii,ee), crystallite_subTstar0_v(:,gg,ii,ee), & + crystallite_Fe, crystallite_Fp, crystallite_Temperature(gg,ii,ee), & + crystallite_subdt(gg,ii,ee), gg, ii, ee) ! collect dot state + enddo; enddo; enddo - stateConverged = crystallite_updateState(g,i,e) ! update state - temperatureConverged = crystallite_updateTemperature(g,i,e) ! update temperature - converged = stateConverged .and. temperatureConverged - endif - if (debugger) then - !$OMP CRITICAL (write2out) - write (6,*) '-------------' - write (6,'(a,x,l,x,l)') 'ontrack + converged:',onTrack,converged - write (6,'(a,/,3(3(f12.4,x)/))') 'pertP of 1 1 1',crystallite_P(1:3,:,g,i,e)/1e6 - write (6,'(a,/,3(3(f12.4,x)/))') 'DP of 1 1 1',(crystallite_P(1:3,:,g,i,e)-myP(1:3,:))/1e6 - write (6,'(a,/,16(6(e12.4,x)/),/,2(f12.4,x))') 'state of 1 1 1',constitutive_state(g,i,e)%p/1e6 - write (6,'(a,/,16(6(e12.4,x)/),/,2(f12.4,x))') 'Dstate of 1 1 1',(constitutive_state(g,i,e)%p-myState)/1e6 - !$OMPEND CRITICAL (write2out) - endif - enddo - if (converged) then ! converged state warrants stiffness update - dPdF_neg(:,:,k,l) = (myP - crystallite_P(:,:,g,i,e))/pert_Fg ! tangent dP_ij/dFg_kl - if (pert_method == 2_pInt) crystallite_dPdF(:,:,k,l,g,i,e) = dPdF_neg(:,:,k,l) + do ee = FEsolving_execElem(1),FEsolving_execElem(2) + myNgrains = homogenization_Ngrains(mesh_element(3,ee)) + do ii = FEsolving_execIP(1,ee),FEsolving_execIP(2,ee) + do gg = 1,myNgrains + if (crystallite_todo(gg,ii,ee)) then + crystallite_stateConverged(gg,ii,ee) = crystallite_updateState(gg,ii,ee) ! update state + crystallite_temperatureConverged(gg,ii,ee) = crystallite_updateTemperature(gg,ii,ee) ! update temperature + crystallite_converged(gg,ii,ee) = crystallite_stateConverged(gg,ii,ee) & + .and. crystallite_temperatureConverged(gg,ii,ee) endif - constitutive_state(g,i,e)%p = myState ! restore unperturbed, converged state, ... - constitutive_dotState(g,i,e)%p = myDotState ! ... dotState, ... - crystallite_Temperature(g,i,e) = myTemperature ! ... temperature, ... - crystallite_Fp(:,:,g,i,e) = myFp ! ... and kinematics - crystallite_invFp(:,:,g,i,e) = myInvFp - crystallite_Fe(:,:,g,i,e) = myFe - crystallite_Lp(:,:,g,i,e) = myLp - crystallite_Tstar_v(:,g,i,e) = myTstar_v - crystallite_P(:,:,g,i,e) = myP - !$OMP CRITICAL (out) - debug_StiffnessStateLoopDistribution(NiterationState) = & - debug_StiffnessstateLoopDistribution(NiterationState) + 1 - !$OMPEND CRITICAL (out) enddo - enddo - endif - if (pert_method == 3_pInt) crystallite_dPdF(:,:,:,:,g,i,e) = 0.5_pReal*(dPdF_neg + dPdF_pos) - else ! grain did not converged - crystallite_dPdF(:,:,:,:,g,i,e) = crystallite_fallbackdPdF(:,:,:,:,g,i,e) ! use (elastic) fallback - endif ! grain convergence - endif ! grain request - enddo ! grain loop - enddo ! ip loop - enddo ! element loop - !$OMPEND PARALLEL DO + enddo + enddo + + 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 ! state loop + + if (all(crystallite_converged(:,:,FEsolving_execELem(1):FEsolving_execElem(2)))) then + crystallite_dPdF(:,:,k,l,g,i,e) = (crystallite_P(:,:,g,i,e) - storedP(:,:,g,i,e))/pert_Fg ! tangent dP_ij/dFg_kl + else ! grain did not converge + crystallite_dPdF(:,:,k,l,g,i,e) = crystallite_fallbackdPdF(:,:,k,l,g,i,e) ! use (elastic) fallback + endif + + do ee = FEsolving_execElem(1),FEsolving_execElem(2) + myNgrains = homogenization_Ngrains(mesh_element(3,ee)) + do ii = FEsolving_execIP(1,ee),FEsolving_execIP(2,ee) + do gg = 1,myNgrains + mySizeState = constitutive_sizeState(gg,ii,ee) + mySizeDotState = constitutive_sizeDotState(gg,ii,ee) + constitutive_state(gg,ii,ee)%p = storedState(1:mySizeState,gg,ii,ee) + constitutive_dotState(gg,ii,ee)%p = storedDotState(1:mySizeDotState,gg,ii,ee) + crystallite_Temperature(gg,ii,ee) = storedTemperature(gg,ii,ee) + crystallite_subF(:,:,gg,ii,ee) = storedF(:,:,gg,ii,ee) + crystallite_Fp(:,:,gg,ii,ee) = storedFp(:,:,gg,ii,ee) + crystallite_invFp(:,:,gg,ii,ee) = storedInvFp(:,:,gg,ii,ee) + crystallite_Fe(:,:,gg,ii,ee) = storedFe(:,:,gg,ii,ee) + crystallite_Lp(:,:,gg,ii,ee) = storedLp(:,:,gg,ii,ee) + crystallite_Tstar_v(:,gg,ii,ee) = storedTstar_v(:,gg,ii,ee) + crystallite_P(:,:,gg,ii,ee) = storedP(:,:,gg,ii,ee) + enddo; enddo; enddo + + enddo; enddo; enddo ! element,ip,grain loop (e,i,g) + crystallite_converged = storedConvergenceFlag + + + + + + endif endif ! jacobian calculation @@ -781,14 +952,14 @@ endsubroutine 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, & - debug_cumDotStateCalls, & - debug_cumDotStateTicks + use debug, only: debugger + use FEsolving, only: cycleCounter, theInc !*** input variables ***! integer(pInt), intent(in):: e, & ! element index @@ -801,21 +972,16 @@ endsubroutine !*** local variables ***! real(pReal), dimension(constitutive_sizeDotState(g,i,e)) :: residuum ! residuum from evolution of microstructure integer(pInt) mySize - integer(pLongInt) tick, & - tock, & - tickrate, & - maxticks + 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 & + + constitutive_previousDotState(g,i,e)%p(1:mySize) * (1.0_pReal-crystallite_statedamper) ! calculate the residuum - call system_clock(count=tick,count_rate=tickrate,count_max=maxticks) - residuum = constitutive_state(g,i,e)%p(1:mySize) - constitutive_subState0(g,i,e)%p(1:mySize) - & - crystallite_subdt(g,i,e) * constitutive_dotState(g,i,e)%p(1:mySize) - call system_clock(count=tock,count_rate=tickrate,count_max=maxticks) - debug_cumDotStateCalls = debug_cumDotStateCalls + 1_pInt - debug_cumDotStateTicks = debug_cumDotStateTicks + tock-tick - if (tock < tick) debug_cumDotStateTicks = debug_cumDotStateTicks + maxticks + 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 NaN occured then return without changing the state if (any(residuum/=residuum)) then @@ -845,9 +1011,13 @@ endsubroutine write(6,*) '::: updateState did not converge',g,i,e endif write(6,*) - write(6,'(a10,/,12(e12.5,x))') 'new state ',constitutive_state(g,i,e)%p(1:mySize) + write(6,'(a,f6.1)') 'crystallite_statedamper',crystallite_statedamper write(6,*) - write(6,'(a,/,12(f12.5,x))') 'resid tolerance',abs(residuum/rTol_crystalliteState/constitutive_state(g,i,e)%p(1:mySize)) + write(6,'(a,/,12(e12.5,x))') 'dotState',constitutive_dotState(g,i,e)%p(1:mySize) + write(6,*) + write(6,'(a,/,12(e12.5,x))') 'new state',constitutive_state(g,i,e)%p(1:mySize) + write(6,*) + write(6,'(a,/,12(f12.1,x))') 'resid tolerance',abs(residuum/rTol_crystalliteState/constitutive_state(g,i,e)%p(1:mySize)) write(6,*) !$OMPEND CRITICAL (write2out) endif @@ -1060,7 +1230,7 @@ LpLoop: do if (NiterationStress > nStress) then if (debugger) then !$OMP CRITICAL (write2out) - write(6,*) '::: integrateStress reached loop limit',g,i,e + write(6,*) '::: integrateStress reached loop limit at ',g,i,e write(6,*) !$OMPEND CRITICAL (write2out) endif @@ -1086,10 +1256,11 @@ LpLoop: do if (tock < tick) debug_cumLpTicks = debug_cumLpTicks + maxticks if (debugger) then !$OMP CRITICAL (write2out) - write(6,*) '::: integrateStress at iteration', NiterationStress + write(6,*) '::: integrateStress at ' ,g,i,e, ' ; iteration ', NiterationStress write(6,*) - write(6,'(a19,3(i3,x),/,3(3(f20.7,x)/))') 'Lp_constitutive at ',g,i,e,Lp_constitutive - write(6,'(a11,3(i3,x),/,3(3(f20.7,x)/))') 'Lpguess at ',g,i,e,Lpguess + write(6,'(a,/,3(3(f20.7,x)/))') 'Lp_constitutive', Lp_constitutive + write(6,'(a,/,3(3(f20.7,x)/))') 'Lpguess', Lpguess + ! call flush(6) !$OMPEND CRITICAL (write2out) endif @@ -1110,7 +1281,7 @@ LpLoop: do if (any(residuum/=residuum) .and. leapfrog == 1.0) then if (debugger) then !$OMP CRITICAL (write2out) - write(6,*) '::: integrateStress encountered NaN at iteration', NiterationStress,'at',g,i,e + write(6,*) '::: integrateStress encountered NaN at ',g,i,e,' ; iteration ', NiterationStress !$OMPEND CRITICAL (write2out) endif return @@ -1144,12 +1315,12 @@ LpLoop: do if (error) then if (debugger) then !$OMP CRITICAL (write2out) - write(6,*) '::: integrateStress failed on dR/dLp inversion at iteration', NiterationStress + write(6,*) '::: integrateStress failed on dR/dLp inversion at ',g,i,e,' ; iteration ', NiterationStress write(6,*) - write(6,'(a9,3(i3,x),/,9(9(f15.3,x)/))') 'dRdLp at ',g,i,e,dRdLp - write(6,'(a20,3(i3,x),/,9(9(f15.3,x)/))') 'dLpdT_constitutive at ',g,i,e,dLpdT_constitutive - write(6,'(a19,3(i3,x),/,3(3(f20.7,x)/))') 'Lp_constitutive at ',g,i,e,Lp_constitutive - write(6,'(a11,3(i3,x),/,3(3(f20.7,x)/))') 'Lpguess at ',g,i,e,Lpguess + write(6,'(a,/,9(9(f15.3,x)/))') 'dRdLp',dRdLp + write(6,'(a,/,9(9(f15.3,x)/))') 'dLpdT_constitutive',dLpdT_constitutive + write(6,'(a,/,3(3(f20.7,x)/))') 'Lp_constitutive',Lp_constitutive + write(6,'(a,/,3(3(f20.7,x)/))') 'Lpguess',Lpguess !$OMPEND CRITICAL (write2out) endif return @@ -1177,7 +1348,7 @@ LpLoop: do if (error) then if (debugger) then !$OMP CRITICAL (write2out) - write(6,*) '::: integrateStress failed on invFp_new inversion at iteration', NiterationStress + write(6,*) '::: integrateStress failed on invFp_new inversion at ',g,i,e,' ; iteration ', NiterationStress write(6,*) write(6,'(a11,3(i3,x),/,3(3(f12.7,x)/))') 'invFp_new at ',g,i,e,invFp_new !$OMPEND CRITICAL (write2out) @@ -1203,7 +1374,7 @@ LpLoop: do crystallite_integrateStress = .true. if (debugger) then !$OMP CRITICAL (write2out) - write(6,*) '::: integrateStress converged at iteration', NiterationStress + write(6,*) '::: integrateStress converged at ',g,i,e,' ; iteration ', NiterationStress write(6,*) write(6,'(a,/,3(3(f12.7,x)/))') 'P / MPa',crystallite_P(:,:,g,i,e)/1e6 write(6,'(a,/,3(3(f12.7,x)/))') 'Lp',crystallite_Lp(:,:,g,i,e) @@ -1257,7 +1428,7 @@ function crystallite_postResults(& !*** local variables ***! real(pReal), dimension(3,3) :: U, R - integer(pInt) k,l,c + integer(pInt) k,l,c logical error c = 0_pInt @@ -1284,8 +1455,9 @@ 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_Temperature(g,i,e), & - dt, crystallite_subdt(g,i,e), g, i, e); c = c+constitutive_sizePostResults(g,i,e) + constitutive_postResults(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, crystallite_Fp, & + crystallite_Temperature(g,i,e), dt, crystallite_subdt(g,i,e), g, i, e) + c = c + constitutive_sizePostResults(g,i,e) return diff --git a/code/debug.f90 b/code/debug.f90 index 159e76200..c2e2acf3d 100644 --- a/code/debug.f90 +++ b/code/debug.f90 @@ -100,7 +100,7 @@ endsubroutine dble(debug_cumLpTicks)*1.0e6_pReal/tickrate/debug_cumLpCalls endif write(6,*) - write(6,'(a33,x,i12)') 'total calls to dotState :',debug_cumDotStateCalls + write(6,'(a33,x,i12)') 'total calls to collectDotState :',debug_cumDotStateCalls if (debug_cumdotStateCalls > 0_pInt) then write(6,'(a33,x,f12.3)') 'total CPU time/s :',dble(debug_cumDotStateTicks)/tickrate write(6,'(a33,x,f12.6)') 'avg CPU time/microsecs per call :',& @@ -191,7 +191,6 @@ endsubroutine enddo write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_MaterialpointLoopDistribution) write(6,*) - call flush(6) endsubroutine