diff --git a/code/CPFEM.f90 b/code/CPFEM.f90 index 0386a7e0c..76fae6c8e 100644 --- a/code/CPFEM.f90 +++ b/code/CPFEM.f90 @@ -73,6 +73,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt use FEsolving, only: FE_init, & parallelExecution, & outdatedFFN1, & + terminallyIll, & cycleCounter, & theInc, & theCycle, & @@ -201,7 +202,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt j = 1:mesh_maxNips, & k = 1:mesh_NcpElems ) & constitutive_state0(i,j,k)%p = constitutive_state(i,j,k)%p ! microstructure of crystallites - write(6,'(a10,/,4(3(f20.8,x),/))') 'aged state',constitutive_state(1,1,1)%p + write(6,'(a10,/,4(3(e20.8,x),/))') 'aged state',constitutive_state(1,1,1)%p do k = 1,mesh_NcpElems do j = 1,mesh_maxNips if (homogenization_sizeState(j,k) > 0_pInt) & @@ -211,10 +212,11 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt endif ! deformation gradient outdated or any actual deformation gradient differs more than relevantStrain from the stored one - if (outdatedFFN1 .or. any(abs(ffn1 - materialpoint_F(:,:,IP,cp_en)) > relevantStrain)) then - if (.not. outdatedFFN1) & + if (terminallyIll .or. outdatedFFN1 .or. any(abs(ffn1 - materialpoint_F(:,:,IP,cp_en)) > relevantStrain)) then + if (.not. terminallyIll .and. .not. outdatedFFN1) then write(6,'(a11,x,i5,x,i2,x,a10,/,3(3(f10.3,x),/))') 'outdated at',cp_en,IP,'FFN1 now:',ffn1(:,1),ffn1(:,2),ffn1(:,3) - outdatedFFN1 = .true. + outdatedFFN1 = .true. + endif CPFEM_cs(1:ngens,IP,cp_en) = CPFEM_odd_stress CPFEM_dcsde(1:ngens,1:ngens,IP,cp_en) = CPFEM_odd_jacobian*math_identity2nd(ngens) @@ -240,25 +242,29 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt CPFEM_calc_done = .true. endif + if (terminallyIll) then + CPFEM_cs(1:ngens,IP,cp_en) = CPFEM_odd_stress + CPFEM_dcsde(1:ngens,1:ngens,IP,cp_en) = CPFEM_odd_jacobian*math_identity2nd(ngens) + else ! translate from P to CS - Kirchhoff = math_mul33x33(materialpoint_P(:,:,IP, cp_en),transpose(materialpoint_F(:,:,IP, cp_en))) - J_inverse = 1.0_pReal/math_det3x3(materialpoint_F(:,:,IP, cp_en)) - CPFEM_cs(1:ngens,IP,cp_en) = math_Mandel33to6(J_inverse*Kirchhoff) + Kirchhoff = math_mul33x33(materialpoint_P(:,:,IP, cp_en),transpose(materialpoint_F(:,:,IP, cp_en))) + J_inverse = 1.0_pReal/math_det3x3(materialpoint_F(:,:,IP, cp_en)) + CPFEM_cs(1:ngens,IP,cp_en) = math_Mandel33to6(J_inverse*Kirchhoff) ! translate from dP/dF to dCS/dE - H = 0.0_pReal - forall(i=1:3,j=1:3,k=1:3,l=1:3,m=1:3,n=1:3) & - H(i,j,k,l) = H(i,j,k,l) + & - materialpoint_F(j,m,IP,cp_en) * & - materialpoint_F(l,n,IP,cp_en) * & - materialpoint_dPdF(i,m,k,n,IP,cp_en) - & - math_I3(j,l)*materialpoint_F(i,m,IP,cp_en)*materialpoint_P(k,m,IP,cp_en) + & - 0.5_pReal*(math_I3(i,k)*Kirchhoff(j,l) + math_I3(j,l)*Kirchhoff(i,k) + & - math_I3(i,l)*Kirchhoff(j,k) + math_I3(j,k)*Kirchhoff(i,l)) - forall(i=1:3,j=1:3,k=1:3,l=1:3) & - H_sym(i,j,k,l)= 0.25_pReal*(H(i,j,k,l)+H(j,i,k,l)+H(i,j,l,k)+H(j,i,l,k)) ! where to use the symmetric version?? - CPFEM_dcsde(1:ngens,1:ngens,IP,cp_en) = math_Mandel3333to66(J_inverse*H) - + H = 0.0_pReal + forall(i=1:3,j=1:3,k=1:3,l=1:3,m=1:3,n=1:3) & + H(i,j,k,l) = H(i,j,k,l) + & + materialpoint_F(j,m,IP,cp_en) * & + materialpoint_F(l,n,IP,cp_en) * & + materialpoint_dPdF(i,m,k,n,IP,cp_en) - & + math_I3(j,l)*materialpoint_F(i,m,IP,cp_en)*materialpoint_P(k,m,IP,cp_en) + & + 0.5_pReal*(math_I3(i,k)*Kirchhoff(j,l) + math_I3(j,l)*Kirchhoff(i,k) + & + math_I3(i,l)*Kirchhoff(j,k) + math_I3(j,k)*Kirchhoff(i,l)) + forall(i=1:3,j=1:3,k=1:3,l=1:3) & + H_sym(i,j,k,l)= 0.25_pReal*(H(i,j,k,l)+H(j,i,k,l)+H(i,j,l,k)+H(j,i,l,k)) ! where to use the symmetric version?? + CPFEM_dcsde(1:ngens,1:ngens,IP,cp_en) = math_Mandel3333to66(J_inverse*H) + endif endif ! --+>> COLLECTION OF FEM DATA AND RETURN OF ODD STRESS AND JACOBIAN <<+-- diff --git a/code/FEsolving.f90 b/code/FEsolving.f90 index ac2a334e6..4a464d530 100644 --- a/code/FEsolving.f90 +++ b/code/FEsolving.f90 @@ -9,7 +9,7 @@ integer(pInt) cycleCounter integer(pInt) theInc,theCycle,theLovl real(pReal) theTime - logical :: lastIncConverged = .false.,outdatedByNewInc = .false.,outdatedFFN1 = .false. + logical :: lastIncConverged = .false.,outdatedByNewInc = .false.,outdatedFFN1 = .false.,terminallyIll = .false. logical :: symmetricSolver = .false. logical :: parallelExecution = .true. integer(pInt), dimension(:,:), allocatable :: FEsolving_execIP diff --git a/code/IO.f90 b/code/IO.f90 index 4ac4cad2b..6dae90389 100644 --- a/code/IO.f90 +++ b/code/IO.f90 @@ -858,7 +858,7 @@ endfunction case (265) msg = 'Limit for crystallite loop too small' case (266) - msg = 'Limit for state loop too small' + msg = 'Limit for crystallite state loop too small' case (267) msg = 'Limit for stress loop too small' case (268) @@ -884,6 +884,9 @@ endfunction case (277) msg = 'Non-positive relevant mismatch in RGC' + case (279) + msg = 'Limit for materialpoint state loop too small' + case (300) msg = 'This material can only be used with elements with three direct stress components' case (500) diff --git a/code/constitutive.f90 b/code/constitutive.f90 index 66b033a88..8b341a745 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -16,11 +16,14 @@ MODULE constitutive type(p_vec), dimension(:,:,:), allocatable :: constitutive_state0, & ! pointer array to microstructure at start of FE inc constitutive_partionedState0, & ! pointer array to microstructure at start of homogenization inc constitutive_subState0, & ! pointer array to microstructure at start of crystallite inc - constitutive_state ! pointer array to current microstructure (end of converged time step) + constitutive_state, & ! pointer array to current microstructure (end of converged time step) + constitutive_dotState ! pointer array to evolution of current microstructure integer(pInt), dimension(:,:,:), allocatable :: constitutive_sizeDotState, & ! size of dotState array constitutive_sizeState, & ! size of state array per grain constitutive_sizePostResults ! size of postResults array per grain - integer(pInt) constitutive_maxSizeDotState,constitutive_maxSizeState,constitutive_maxSizePostResults + integer(pInt) constitutive_maxSizeDotState, & + constitutive_maxSizeState, & + constitutive_maxSizePostResults CONTAINS !**************************************** @@ -28,8 +31,8 @@ CONTAINS !* - constitutive_homogenizedC !* - constitutive_microstructure !* - constitutive_LpAndItsTangent -!* - constitutive_dotState -!* - constitutive_dotTemperature +!* - constitutive_collectDotState +!* - constitutive_collectDotTemperature !* - constitutive_postResults !**************************************** @@ -46,6 +49,7 @@ subroutine constitutive_init() use constitutive_j2 use constitutive_phenopowerlaw use constitutive_dislobased + use constitutive_nonlocal integer(pInt), parameter :: fileunit = 200 integer(pInt) e,i,g,p,myInstance,myNgrains @@ -58,6 +62,7 @@ subroutine constitutive_init() call constitutive_j2_init(fileunit) ! parse all phases of this constitution call constitutive_phenopowerlaw_init(fileunit) call constitutive_dislobased_init(fileunit) + call constitutive_nonlocal_init(fileunit) close(fileunit) @@ -78,6 +83,9 @@ subroutine constitutive_init() case (constitutive_dislobased_label) thisOutput => constitutive_dislobased_output thisSize => constitutive_dislobased_sizePostResult + case (constitutive_nonlocal_label) + thisOutput => constitutive_nonlocal_output + thisSize => constitutive_nonlocal_sizePostResult case default knownConstitution = .false. end select @@ -101,10 +109,10 @@ subroutine constitutive_init() allocate(constitutive_partionedState0(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) allocate(constitutive_subState0(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) allocate(constitutive_state(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) - allocate(constitutive_sizeDotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_sizeDotState = 0_pInt - allocate(constitutive_sizeState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_sizeState = 0_pInt - allocate(constitutive_sizePostResults(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_sizePostResults = 0_pInt - + allocate(constitutive_dotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) + allocate(constitutive_sizeDotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_sizeDotState = 0_pInt + allocate(constitutive_sizeState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_sizeState = 0_pInt + allocate(constitutive_sizePostResults(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)); constitutive_sizePostResults = 0_pInt do e = 1,mesh_NcpElems ! loop over elements myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = 1,FE_Nips(mesh_element(2,e)) ! loop over IPs @@ -117,28 +125,41 @@ subroutine constitutive_init() allocate(constitutive_partionedState0(g,i,e)%p(constitutive_j2_sizeState(myInstance))) allocate(constitutive_subState0(g,i,e)%p(constitutive_j2_sizeState(myInstance))) allocate(constitutive_state(g,i,e)%p(constitutive_j2_sizeState(myInstance))) + allocate(constitutive_dotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) constitutive_state0(g,i,e)%p = constitutive_j2_stateInit(myInstance) - constitutive_sizeDotState(g,i,e) = constitutive_j2_sizeDotState(myInstance) constitutive_sizeState(g,i,e) = constitutive_j2_sizeState(myInstance) + constitutive_sizeDotState(g,i,e) = constitutive_j2_sizeDotState(myInstance) constitutive_sizePostResults(g,i,e) = constitutive_j2_sizePostResults(myInstance) - case (constitutive_phenopowerlaw_label) + case (constitutive_phenopowerlaw_label) allocate(constitutive_state0(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) allocate(constitutive_partionedState0(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) allocate(constitutive_subState0(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) allocate(constitutive_state(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) + allocate(constitutive_dotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) constitutive_state0(g,i,e)%p = constitutive_phenopowerlaw_stateInit(myInstance) - constitutive_sizeDotState(g,i,e) = constitutive_phenopowerlaw_sizeDotState(myInstance) constitutive_sizeState(g,i,e) = constitutive_phenopowerlaw_sizeState(myInstance) + constitutive_sizeDotState(g,i,e) = constitutive_phenopowerlaw_sizeDotState(myInstance) constitutive_sizePostResults(g,i,e) = constitutive_phenopowerlaw_sizePostResults(myInstance) case (constitutive_dislobased_label) allocate(constitutive_state0(g,i,e)%p(constitutive_dislobased_sizeState(myInstance))) allocate(constitutive_partionedState0(g,i,e)%p(constitutive_dislobased_sizeState(myInstance))) allocate(constitutive_subState0(g,i,e)%p(constitutive_dislobased_sizeState(myInstance))) allocate(constitutive_state(g,i,e)%p(constitutive_dislobased_sizeState(myInstance))) + allocate(constitutive_dotState(g,i,e)%p(constitutive_dislobased_sizeDotState(myInstance))) constitutive_state0(g,i,e)%p = constitutive_dislobased_stateInit(myInstance) - constitutive_sizeDotState(g,i,e) = constitutive_dislobased_sizeDotState(myInstance) constitutive_sizeState(g,i,e) = constitutive_dislobased_sizeState(myInstance) + constitutive_sizeDotState(g,i,e) = constitutive_dislobased_sizeDotState(myInstance) constitutive_sizePostResults(g,i,e) = constitutive_dislobased_sizePostResults(myInstance) + case (constitutive_nonlocal_label) + allocate(constitutive_state0(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance))) + allocate(constitutive_partionedState0(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance))) + allocate(constitutive_subState0(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance))) + allocate(constitutive_state(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance))) + allocate(constitutive_dotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) + constitutive_state0(g,i,e)%p = constitutive_nonlocal_stateInit(myInstance) + constitutive_sizeState(g,i,e) = constitutive_nonlocal_sizeState(myInstance) + constitutive_sizeDotState(g,i,e) = constitutive_nonlocal_sizeDotState(myInstance) + constitutive_sizePostResults(g,i,e) = constitutive_nonlocal_sizePostResults(myInstance) case default call IO_error(200,material_phase(g,i,e)) ! unknown constitution end select @@ -146,8 +167,9 @@ subroutine constitutive_init() enddo enddo enddo - constitutive_maxSizeDotState = maxval(constitutive_sizeDotState) + constitutive_maxSizeState = maxval(constitutive_sizeState) + constitutive_maxSizeDotState = maxval(constitutive_sizeDotState) constitutive_maxSizePostResults = maxval(constitutive_sizePostResults) write(6,*) @@ -157,6 +179,7 @@ subroutine constitutive_init() write(6,'(a32,x,7(i5,x))') 'constitutive_partionedState0: ', shape(constitutive_partionedState0) write(6,'(a32,x,7(i5,x))') 'constitutive_subState0: ', shape(constitutive_subState0) write(6,'(a32,x,7(i5,x))') 'constitutive_state: ', shape(constitutive_state) + write(6,'(a32,x,7(i5,x))') 'constitutive_dotState: ', shape(constitutive_dotState) 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) @@ -166,7 +189,7 @@ subroutine constitutive_init() return -end subroutine +endsubroutine function constitutive_homogenizedC(ipc,ip,el) @@ -183,6 +206,7 @@ function constitutive_homogenizedC(ipc,ip,el) use constitutive_j2 use constitutive_phenopowerlaw use constitutive_dislobased + use constitutive_nonlocal implicit none !* Definition of variables @@ -190,19 +214,26 @@ function constitutive_homogenizedC(ipc,ip,el) real(pReal), dimension(6,6) :: constitutive_homogenizedC select case (phase_constitution(material_phase(ipc,ip,el))) + case (constitutive_j2_label) constitutive_homogenizedC = constitutive_j2_homogenizedC(constitutive_state,ipc,ip,el) + case (constitutive_phenopowerlaw_label) constitutive_homogenizedC = constitutive_phenopowerlaw_homogenizedC(constitutive_state,ipc,ip,el) + case (constitutive_dislobased_label) constitutive_homogenizedC = constitutive_dislobased_homogenizedC(constitutive_state,ipc,ip,el) + + case (constitutive_nonlocal_label) + constitutive_homogenizedC = constitutive_nonlocal_homogenizedC(constitutive_state,ipc,ip,el) + end select return -end function +endfunction -subroutine constitutive_microstructure(Temperature,ipc,ip,el) +subroutine constitutive_microstructure(Temperature,Fp,ipc,ip,el) !********************************************************************* !* This function calculates from state needed variables * !* INPUT: * @@ -212,30 +243,43 @@ subroutine constitutive_microstructure(Temperature,ipc,ip,el) !* - 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 material, only: phase_constitution, & + material_phase, & + homogenization_maxNgrains + use mesh, only: mesh_NcpElems, & + mesh_maxNips use constitutive_j2 use constitutive_phenopowerlaw use constitutive_dislobased + use constitutive_nonlocal implicit none !* Definition of variables -integer(pInt) ipc,ip,el -real(pReal) Temperature +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) :: Fp select case (phase_constitution(material_phase(ipc,ip,el))) + case (constitutive_j2_label) call constitutive_j2_microstructure(Temperature,constitutive_state,ipc,ip,el) + case (constitutive_phenopowerlaw_label) call constitutive_phenopowerlaw_microstructure(Temperature,constitutive_state,ipc,ip,el) + case (constitutive_dislobased_label) call constitutive_dislobased_microstructure(Temperature,constitutive_state,ipc,ip,el) + + case (constitutive_nonlocal_label) + call constitutive_nonlocal_microstructure(Temperature, Fp, constitutive_state,ipc,ip,el) + end select -end subroutine +endsubroutine -subroutine constitutive_LpAndItsTangent(Lp,dLp_dTstar, Tstar_v,Temperature,ipc,ip,el) +subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, Temperature, ipc, ip, el) !********************************************************************* !* This subroutine contains the constitutive equation for * !* calculating the velocity gradient * @@ -253,6 +297,7 @@ subroutine constitutive_LpAndItsTangent(Lp,dLp_dTstar, Tstar_v,Temperature,ipc,i use constitutive_j2 use constitutive_phenopowerlaw use constitutive_dislobased + use constitutive_nonlocal implicit none !* Definition of variables @@ -263,19 +308,26 @@ subroutine constitutive_LpAndItsTangent(Lp,dLp_dTstar, Tstar_v,Temperature,ipc,i real(pReal), dimension(9,9) :: dLp_dTstar select case (phase_constitution(material_phase(ipc,ip,el))) + case (constitutive_j2_label) call constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el) + case (constitutive_phenopowerlaw_label) call constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el) + case (constitutive_dislobased_label) call constitutive_dislobased_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el) + + case (constitutive_nonlocal_label) + call constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, Temperature, constitutive_state, ipc, ip, el) + end select return -end subroutine +endsubroutine -function constitutive_dotState(Tstar_v,Temperature,ipc,ip,el) +subroutine constitutive_collectDotState(Tstar_v, Fp, invFp, Temperature, ipc, ip, el) !********************************************************************* !* This subroutine contains the constitutive equation for * !* calculating the rate of change of microstructure * @@ -289,28 +341,37 @@ function constitutive_dotState(Tstar_v,Temperature,ipc,ip,el) !* - constitutive_dotState : evolution of state variable * !********************************************************************* use prec, only: pReal,pInt + use debug use material, only: phase_constitution,material_phase use constitutive_j2 use constitutive_phenopowerlaw use constitutive_dislobased + use constitutive_nonlocal implicit none !* Definition of variables integer(pInt) ipc,ip,el real(pReal) Temperature + real(pReal), dimension(3,3) :: Fp, invFp real(pReal), dimension(6) :: Tstar_v - real(pReal), dimension(constitutive_sizeDotState(ipc,ip,el)) :: constitutive_dotState select case (phase_constitution(material_phase(ipc,ip,el))) + case (constitutive_j2_label) - constitutive_dotState = constitutive_j2_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el) + constitutive_dotState(ipc,ip,el)%p = constitutive_j2_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el) + case (constitutive_phenopowerlaw_label) - constitutive_dotState = constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el) + constitutive_dotState(ipc,ip,el)%p = constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el) + case (constitutive_dislobased_label) - constitutive_dotState = constitutive_dislobased_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el) + constitutive_dotState(ipc,ip,el)%p = constitutive_dislobased_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el) + + case (constitutive_nonlocal_label) + call constitutive_nonlocal_dotState(constitutive_dotState,Tstar_v,Fp,invFp,Temperature,constitutive_state,ipc,ip,el) + end select return -end function +endsubroutine function constitutive_dotTemperature(Tstar_v,Temperature,ipc,ip,el) @@ -331,25 +392,32 @@ function constitutive_dotTemperature(Tstar_v,Temperature,ipc,ip,el) use constitutive_j2 use constitutive_phenopowerlaw use constitutive_dislobased + use constitutive_nonlocal implicit none !* Definition of variables integer(pInt) ipc,ip,el real(pReal) Temperature - real(pReal), dimension(6) :: Tstar_v real(pReal) constitutive_dotTemperature + real(pReal), dimension(6) :: Tstar_v select case (phase_constitution(material_phase(ipc,ip,el))) + case (constitutive_j2_label) constitutive_dotTemperature = constitutive_j2_dotTemperature(Tstar_v,Temperature,constitutive_state,ipc,ip,el) + case (constitutive_phenopowerlaw_label) constitutive_dotTemperature = constitutive_phenopowerlaw_dotTemperature(Tstar_v,Temperature,constitutive_state,ipc,ip,el) + case (constitutive_dislobased_label) constitutive_dotTemperature = constitutive_dislobased_dotTemperature(Tstar_v,Temperature,constitutive_state,ipc,ip,el) - + + case (constitutive_nonlocal_label) + constitutive_dotTemperature = constitutive_nonlocal_dotTemperature(Tstar_v,Temperature,constitutive_state,ipc,ip,el) + end select return -end function +endfunction pure function constitutive_postResults(Tstar_v,Temperature,dt,ipc,ip,el) @@ -367,6 +435,7 @@ pure function constitutive_postResults(Tstar_v,Temperature,dt,ipc,ip,el) use constitutive_j2 use constitutive_phenopowerlaw use constitutive_dislobased + use constitutive_nonlocal implicit none !* Definition of variables @@ -377,17 +446,23 @@ pure function constitutive_postResults(Tstar_v,Temperature,dt,ipc,ip,el) constitutive_postResults = 0.0_pReal select case (phase_constitution(material_phase(ipc,ip,el))) + case (constitutive_j2_label) constitutive_postResults = constitutive_j2_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el) + case (constitutive_phenopowerlaw_label) constitutive_postResults = constitutive_phenopowerlaw_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el) + case (constitutive_dislobased_label) constitutive_postResults = constitutive_dislobased_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el) - + + case (constitutive_nonlocal_label) + constitutive_postResults = constitutive_nonlocal_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el) + end select return -end function +endfunction END MODULE \ No newline at end of file diff --git a/code/constitutive_nonlocal.f90 b/code/constitutive_nonlocal.f90 new file mode 100644 index 000000000..8bbda8251 --- /dev/null +++ b/code/constitutive_nonlocal.f90 @@ -0,0 +1,1113 @@ + +!************************************ +!* Module: CONSTITUTIVE_NONLOCAL * +!************************************ +!* contains: * +!* - constitutive equations * +!* - parameters definition * +!************************************ + + +MODULE constitutive_nonlocal + +!* Include other modules +use prec, only: pReal,pInt +implicit none + + +!* Definition of parameters +character (len=*), parameter :: constitutive_nonlocal_label = 'nonlocal' +character(len=16), dimension(7), parameter :: constitutive_nonlocal_stateList = (/ 'rhoEdgePos ', & + 'rhoEdgeNeg ', & + 'rhoScrewPos ', & + 'rhoScrewNeg ', & + 'rhoForest ', & + 'tauSlipThreshold', & + 'backStress_v ' /) ! list of microstructural state variables +character(len=16), dimension(4), parameter :: constitutive_nonlocal_stateListBasic = constitutive_nonlocal_stateList(1:4) ! list of "basic" microstructural state variables that are independent from other state variables +character(len=16), dimension(3), parameter :: constitutive_nonlocal_stateListDependent = constitutive_nonlocal_stateList(5:7) ! list of microstructural state variables that depend on other state variables +real(pReal), parameter :: kB = 1.38e-23_pReal ! Physical parameter, Boltzmann constant in J/Kelvin + + +!* Definition of global variables +integer(pInt), dimension(:), allocatable :: constitutive_nonlocal_sizeDotState, & ! number of dotStates + constitutive_nonlocal_sizeState, & ! total number of microstructural state variables + constitutive_nonlocal_sizePostResults ! cumulative size of post results +integer(pInt), dimension(:,:), allocatable, target :: constitutive_nonlocal_sizePostResult ! size of each post result output +character(len=64), dimension(:,:), allocatable, target :: constitutive_nonlocal_output ! name of each post result output + +character(len=32), dimension(:), allocatable :: constitutive_nonlocal_structureName ! name of the lattice structure +integer(pInt), dimension(:), allocatable :: constitutive_nonlocal_structure, & ! number representing the kind of lattice structure + constitutive_nonlocal_totalNslip ! total number of active slip systems for each instance +integer(pInt), dimension(:,:), allocatable :: constitutive_nonlocal_Nslip, & ! number of active slip systems for each family and instance + constitutive_nonlocal_slipFamily, & ! lookup table relating active slip system to slip family for each instance + constitutive_nonlocal_slipSystemLattice ! lookup table relating active slip system index to lattice slip system index for each instance + +real(pReal), dimension(:), allocatable :: constitutive_nonlocal_CoverA, & ! c/a ratio for hex type lattice + constitutive_nonlocal_C11, & ! C11 element in elasticity matrix + constitutive_nonlocal_C12, & ! C12 element in elasticity matrix + constitutive_nonlocal_C13, & ! C13 element in elasticity matrix + constitutive_nonlocal_C33, & ! C33 element in elasticity matrix + constitutive_nonlocal_C44, & ! C44 element in elasticity matrix + constitutive_nonlocal_Gmod, & ! shear modulus + constitutive_nonlocal_nu ! poisson's ratio +real(pReal), dimension(:,:,:), allocatable :: constitutive_nonlocal_Cslip_66 ! elasticity matrix in Mandel notation for each instance +real(pReal), dimension(:,:,:,:,:), allocatable :: constitutive_nonlocal_Cslip_3333 ! elasticity matrix for each instance +real(pReal), dimension(:,:), allocatable :: constitutive_nonlocal_rhoEdgePos0, & ! initial edge_pos dislocation density per slip system for each family and instance + constitutive_nonlocal_rhoEdgeNeg0, & ! initial edge_neg dislocation density per slip system for each family and instance + constitutive_nonlocal_rhoScrewPos0, & ! initial screw_pos dislocation density per slip system for each family and instance + constitutive_nonlocal_rhoScrewNeg0, & ! initial screw_neg dislocation density per slip system for each family and instance + constitutive_nonlocal_v0BySlipFamily, & ! initial dislocation velocity [m/s] for each family and instance + constitutive_nonlocal_v0BySlipSystem, & ! initial dislocation velocity [m/s] for each slip system and instance + constitutive_nonlocal_burgersBySlipFamily, & ! absolute length of burgers vector [m] for each family and instance + constitutive_nonlocal_burgersBySlipSystem, & ! absolute length of burgers vector [m] for each slip system and instance + constitutive_nonlocal_interactionSlipSlip ! coefficients for slip-slip interaction for each interaction type and instance + +real(pReal), dimension(:,:,:), allocatable :: constitutive_nonlocal_forestProjectionEdge, & ! matrix of forest projections of edge dislocations for each instance + 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 + + +CONTAINS +!**************************************** +!* - constitutive_init +!* - constitutive_stateInit +!* - constitutive_homogenizedC +!* - constitutive_microstructure +!* - constitutive_LpAndItsTangent +!* - constitutive_dotState +!* - constitutive_dotTemperature +!* - constitutive_postResults +!**************************************** + + +!************************************** +!* Module initialization * +!************************************** +subroutine constitutive_nonlocal_init(file) + +use prec, only: pInt, pReal +use math, only: math_Mandel3333to66, & + math_Voigt66to3333, & + math_mul3x3 +use IO, only: IO_lc, & + IO_getTag, & + IO_isBlank, & + IO_stringPos, & + IO_stringValue, & + IO_floatValue, & + IO_intValue, & + IO_error +use material, only: phase_constitution, & + phase_constitutionInstance, & + phase_Noutput +use lattice, only: lattice_maxNslipFamily, & + lattice_maxNtwinFamily, & + lattice_maxNslip, & + lattice_maxNtwin, & + lattice_maxNinteraction, & + lattice_NslipSystem, & + lattice_NtwinSystem, & + lattice_initializeStructure, & + lattice_Qtwin, & + lattice_sd, & + lattice_sn, & + lattice_st, & + lattice_interactionSlipSlip + +!*** output variables + +!*** input variables +integer(pInt), intent(in) :: file + +!*** local variables +integer(pInt), parameter :: maxNchunks = 21 +integer(pInt), dimension(1+2*maxNchunks) :: positions +integer(pInt) section, & + maxNinstance, & + maxTotalNslip, & + myStructure, & + f, & ! index of my slip family + i, & ! index of my instance of this constitution + j, & + k, & + l, & + o, & ! index of my output + s, & ! index of my slip system + s1, & ! index of my slip system + s2, & ! index of my slip system + output, & + mySize +character(len=64) tag +character(len=1024) line + + +write(6,*) +write(6,'(a20,a20,a12)') '<<<+- constitutive_',constitutive_nonlocal_label,' init -+>>>' +write(6,*) + +maxNinstance = count(phase_constitution == constitutive_nonlocal_label) +if (maxNinstance == 0) return ! we don't have to do anything if there's no instance for this constitutive law + + +!*** space allocation for global variables + +allocate(constitutive_nonlocal_sizeDotState(maxNinstance)); constitutive_nonlocal_sizeDotState = 0_pInt +allocate(constitutive_nonlocal_sizeState(maxNinstance)); constitutive_nonlocal_sizeState = 0_pInt +allocate(constitutive_nonlocal_sizePostResults(maxNinstance)); constitutive_nonlocal_sizePostResults = 0_pInt +allocate(constitutive_nonlocal_sizePostResult(maxval(phase_Noutput), maxNinstance));constitutive_nonlocal_sizePostResult = 0_pInt +allocate(constitutive_nonlocal_output(maxval(phase_Noutput), maxNinstance)); constitutive_nonlocal_output = '' + +allocate(constitutive_nonlocal_structureName(maxNinstance)); constitutive_nonlocal_structureName = '' +allocate(constitutive_nonlocal_structure(maxNinstance)); constitutive_nonlocal_structure = 0_pInt +allocate(constitutive_nonlocal_Nslip(lattice_maxNslipFamily, maxNinstance)); constitutive_nonlocal_Nslip = 0_pInt +allocate(constitutive_nonlocal_slipFamily(lattice_maxNslip, maxNinstance)); constitutive_nonlocal_slipFamily = 0_pInt +allocate(constitutive_nonlocal_slipSystemLattice(lattice_maxNslip, maxNinstance)); constitutive_nonlocal_slipSystemLattice = 0_pInt +allocate(constitutive_nonlocal_totalNslip(maxNinstance)); constitutive_nonlocal_totalNslip = 0_pInt + +allocate(constitutive_nonlocal_CoverA(maxNinstance)); constitutive_nonlocal_CoverA = 0.0_pReal +allocate(constitutive_nonlocal_C11(maxNinstance)); constitutive_nonlocal_C11 = 0.0_pReal +allocate(constitutive_nonlocal_C12(maxNinstance)); constitutive_nonlocal_C12 = 0.0_pReal +allocate(constitutive_nonlocal_C13(maxNinstance)); constitutive_nonlocal_C13 = 0.0_pReal +allocate(constitutive_nonlocal_C33(maxNinstance)); constitutive_nonlocal_C33 = 0.0_pReal +allocate(constitutive_nonlocal_C44(maxNinstance)); constitutive_nonlocal_C44 = 0.0_pReal +allocate(constitutive_nonlocal_Gmod(maxNinstance)); constitutive_nonlocal_Gmod = 0.0_pReal +allocate(constitutive_nonlocal_nu(maxNinstance)); constitutive_nonlocal_nu = 0.0_pReal +allocate(constitutive_nonlocal_Cslip_66(6,6,maxNinstance)); constitutive_nonlocal_Cslip_66 = 0.0_pReal +allocate(constitutive_nonlocal_Cslip_3333(3,3,3,3,maxNinstance)); constitutive_nonlocal_Cslip_3333 = 0.0_pReal + +allocate(constitutive_nonlocal_rhoEdgePos0(lattice_maxNslipFamily, maxNinstance)); constitutive_nonlocal_rhoEdgePos0 = 0.0_pReal +allocate(constitutive_nonlocal_rhoEdgeNeg0(lattice_maxNslipFamily, maxNinstance)); constitutive_nonlocal_rhoEdgeNeg0 = 0.0_pReal +allocate(constitutive_nonlocal_rhoScrewPos0(lattice_maxNslipFamily, maxNinstance)); constitutive_nonlocal_rhoScrewPos0 = 0.0_pReal +allocate(constitutive_nonlocal_rhoScrewNeg0(lattice_maxNslipFamily, maxNinstance)); constitutive_nonlocal_rhoScrewNeg0 = 0.0_pReal +allocate(constitutive_nonlocal_v0BySlipFamily(lattice_maxNslipFamily, maxNinstance)); + constitutive_nonlocal_v0BySlipFamily = 0.0_pReal +allocate(constitutive_nonlocal_burgersBySlipFamily(lattice_maxNslipFamily, maxNinstance)); + constitutive_nonlocal_burgersBySlipFamily = 0.0_pReal + +allocate(constitutive_nonlocal_interactionSlipSlip(lattice_maxNinteraction, maxNinstance)) + constitutive_nonlocal_interactionSlipSlip = 0.0_pReal + + +!*** readout data from material.config file + +rewind(file) +line = '' +section = 0 + +do while (IO_lc(IO_getTag(line,'<','>')) /= 'phase') ! wind forward to + read(file,'(a1024)',END=100) line +enddo + +do ! read thru sections of phase part + read(file,'(a1024)',END=100) line + if (IO_isBlank(line)) cycle ! skip empty lines + if (IO_getTag(line,'<','>') /= '') exit ! stop at next part + if (IO_getTag(line,'[',']') /= '') then ! next section + section = section + 1 + output = 0 ! reset output counter + endif + if (section > 0 .and. phase_constitution(section) == constitutive_nonlocal_label) then ! one of my sections + i = phase_constitutionInstance(section) ! which instance of my constitution is present phase + positions = IO_stringPos(line,maxNchunks) + tag = IO_lc(IO_stringValue(line,positions,1)) ! extract key + select case(tag) + case ('(output)') + output = output + 1 + constitutive_nonlocal_output(output,i) = IO_lc(IO_stringValue(line,positions,2)) + case ('lattice_structure') + constitutive_nonlocal_structureName(i) = IO_lc(IO_stringValue(line,positions,2)) + case ('covera_ratio') + constitutive_nonlocal_CoverA(i) = IO_floatValue(line,positions,2) + case ('c11') + constitutive_nonlocal_C11(i) = IO_floatValue(line,positions,2) + case ('c12') + constitutive_nonlocal_C12(i) = IO_floatValue(line,positions,2) + case ('c13') + constitutive_nonlocal_C13(i) = IO_floatValue(line,positions,2) + case ('c33') + constitutive_nonlocal_C33(i) = IO_floatValue(line,positions,2) + case ('c44') + constitutive_nonlocal_C44(i) = IO_floatValue(line,positions,2) + case ('nslip') + forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_Nslip(f,i) = IO_intValue(line,positions,1+f) + case ('rhoedgepos0') + forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_rhoEdgePos0(f,i) = IO_floatValue(line,positions,1+f) + case ('rhoedgeneg0') + forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_rhoEdgeNeg0(f,i) = IO_floatValue(line,positions,1+f) + case ('rhoscrewpos0') + forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_rhoScrewPos0(f,i) = IO_floatValue(line,positions,1+f) + case ('rhoscrewneg0') + forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_rhoScrewNeg0(f,i) = IO_floatValue(line,positions,1+f) + case ('v0') + forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_v0BySlipFamily(f,i) = IO_floatValue(line,positions,1+f) + case ('burgers') + forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_burgersBySlipFamily(f,i) = IO_floatValue(line,positions,1+f) + case ('interaction_slipslip') + forall (f = 1:lattice_maxNinteraction) constitutive_nonlocal_interactionSlipSlip(f,i) = IO_floatValue(line,positions,1+f) + end select + endif +enddo + + +100 do i = 1,maxNinstance + + constitutive_nonlocal_structure(i) = & + lattice_initializeStructure(constitutive_nonlocal_structureName(i), constitutive_nonlocal_CoverA(i)) ! our lattice structure is defined in the material.config file by the structureName (and the c/a ratio) + + +!*** sanity checks + + if (constitutive_nonlocal_structure(i) < 1 .or. constitutive_nonlocal_structure(i) > 3) call IO_error(205) + + if (sum(constitutive_nonlocal_Nslip(:,i)) <= 0) call IO_error(225) + + if (any( constitutive_nonlocal_rhoEdgePos0(:,i) < 0.0_pReal & + .and. constitutive_nonlocal_Nslip(:,i) > 0)) call IO_error(220) + + if (any( constitutive_nonlocal_rhoEdgeNeg0(:,i) < 0.0_pReal & + .and. constitutive_nonlocal_Nslip(:,i) > 0)) call IO_error(220) + + if (any( constitutive_nonlocal_rhoScrewPos0(:,i) < 0.0_pReal & + .and. constitutive_nonlocal_Nslip(:,i) > 0)) call IO_error(220) + + if (any( constitutive_nonlocal_rhoScrewNeg0(:,i) < 0.0_pReal & + .and. constitutive_nonlocal_Nslip(:,i) > 0)) call IO_error(220) + + if (any( constitutive_nonlocal_burgersBySlipFamily(:,i) <= 0.0_pReal & + .and. constitutive_nonlocal_Nslip(:,i) > 0)) call IO_error(221) + + if (any( constitutive_nonlocal_v0BySlipFamily(:,i) <= 0.0_pReal & + .and. constitutive_nonlocal_Nslip(:,i) > 0)) call IO_error(-1) + + +!*** determine total number of active slip systems + + constitutive_nonlocal_Nslip(:,i) & + = min( lattice_NslipSystem(:, constitutive_nonlocal_structure(i)), constitutive_nonlocal_Nslip(:,i) ) ! we can't use more slip systems per family than specified in lattice + constitutive_nonlocal_totalNslip(i) = sum(constitutive_nonlocal_Nslip(:,i)) + +enddo + + +!*** allocation of variables whose size depends on the total number of active slip systems + +maxTotalNslip = maxval(constitutive_nonlocal_totalNslip) +allocate(constitutive_nonlocal_burgersBySlipSystem(maxTotalNslip, maxNinstance)) + constitutive_nonlocal_burgersBySlipSystem = 0.0_pReal +allocate(constitutive_nonlocal_v0BySlipSystem(maxTotalNslip, maxNinstance)) + constitutive_nonlocal_v0BySlipSystem = 0.0_pReal +allocate(constitutive_nonlocal_forestProjectionEdge(maxTotalNslip, maxTotalNslip, maxNinstance)) + constitutive_nonlocal_forestProjectionEdge = 0.0_pReal +allocate(constitutive_nonlocal_forestProjectionScrew(maxTotalNslip, maxTotalNslip, maxNinstance)) + constitutive_nonlocal_forestProjectionScrew = 0.0_pReal +allocate(constitutive_nonlocal_interactionMatrixSlipSlip(maxTotalNslip, maxTotalNslip, maxNinstance)) + constitutive_nonlocal_interactionMatrixSlipSlip = 0.0_pReal + +do i = 1,maxNinstance + + myStructure = constitutive_nonlocal_structure(i) ! lattice structure of this instance + +!*** Inverse lookup of my slip system family and the slip system in lattice + + l = 0_pInt + do f = 1,lattice_maxNslipFamily + do s = 1,constitutive_nonlocal_Nslip(f,i) + l = l + 1 + constitutive_nonlocal_slipFamily(l,i) = f + constitutive_nonlocal_slipSystemLattice(l,i) = sum(lattice_NslipSystem(1:f-1, myStructure)) + s + enddo; enddo + + +!*** determine size of state array + + constitutive_nonlocal_sizeState(i) = (size(constitutive_nonlocal_stateList)-1) * constitutive_nonlocal_totalNslip(i) + 6_pInt ! the size of the list of states times the number of active slip systems gives the required size for the state array + constitutive_nonlocal_sizeDotState(i) = size(constitutive_nonlocal_stateListBasic) * constitutive_nonlocal_totalNslip(i) ! the size of the list of basic states times the number of active slip systems gives the required size for the dotState array + + +!*** determine size of postResults array + + do o = 1,maxval(phase_Noutput) + select case(constitutive_nonlocal_output(o,i)) + case( 'dislocationdensity', & + 'shearrate_slip', & + 'resolvedstress_slip', & + 'resistance_slip') + mySize = constitutive_nonlocal_totalNslip(i) + case default + mySize = 0_pInt + end select + + if (mySize > 0_pInt) then ! any meaningful output found + constitutive_nonlocal_sizePostResult(o,i) = mySize + constitutive_nonlocal_sizePostResults(i) = constitutive_nonlocal_sizePostResults(i) + mySize + endif + enddo + + +!*** elasticity matrix and shear modulus according to material.config + + select case (myStructure) + case(1:2) ! cubic(s) + forall(k=1:3) + forall(j=1:3) constitutive_nonlocal_Cslip_66(k,j,i) = constitutive_nonlocal_C12(i) + constitutive_nonlocal_Cslip_66(k,k,i) = constitutive_nonlocal_C11(i) + constitutive_nonlocal_Cslip_66(k+3,k+3,i) = constitutive_nonlocal_C44(i) + end forall + case(3:) ! all hex + constitutive_nonlocal_Cslip_66(1,1,i) = constitutive_nonlocal_C11(i) + constitutive_nonlocal_Cslip_66(2,2,i) = constitutive_nonlocal_C11(i) + constitutive_nonlocal_Cslip_66(3,3,i) = constitutive_nonlocal_C33(i) + constitutive_nonlocal_Cslip_66(1,2,i) = constitutive_nonlocal_C12(i) + constitutive_nonlocal_Cslip_66(2,1,i) = constitutive_nonlocal_C12(i) + constitutive_nonlocal_Cslip_66(1,3,i) = constitutive_nonlocal_C13(i) + constitutive_nonlocal_Cslip_66(3,1,i) = constitutive_nonlocal_C13(i) + constitutive_nonlocal_Cslip_66(2,3,i) = constitutive_nonlocal_C13(i) + constitutive_nonlocal_Cslip_66(3,2,i) = constitutive_nonlocal_C13(i) + constitutive_nonlocal_Cslip_66(4,4,i) = constitutive_nonlocal_C44(i) + constitutive_nonlocal_Cslip_66(5,5,i) = constitutive_nonlocal_C44(i) + constitutive_nonlocal_Cslip_66(6,6,i) = 0.5_pReal*(constitutive_nonlocal_C11(i)- constitutive_nonlocal_C12(i)) + end select + constitutive_nonlocal_Cslip_66(:,:,i) = math_Mandel3333to66(math_Voigt66to3333(constitutive_nonlocal_Cslip_66(:,:,i))) + constitutive_nonlocal_Cslip_3333(:,:,:,:,i) = math_Voigt66to3333(constitutive_nonlocal_Cslip_66(:,:,i)) + + constitutive_nonlocal_Gmod(i) = constitutive_nonlocal_C44(i) ! shear modulus is given by elastic constant C44 + constitutive_nonlocal_nu(i) = constitutive_nonlocal_C12(i) / constitutive_nonlocal_C11(i) + + +!*** burgers vector and initial dislocation velocity for each slip system + + do s = 1,constitutive_nonlocal_totalNslip(i) + constitutive_nonlocal_burgersBySlipSystem(s,i) & + = constitutive_nonlocal_burgersBySlipFamily( constitutive_nonlocal_slipFamily(s,i), i ) + constitutive_nonlocal_v0BySlipSystem(s,i) = constitutive_nonlocal_v0BySlipFamily(constitutive_nonlocal_slipFamily(s,i),i) + enddo + + +!*** calculation of forest projections for edge and screw dislocations + + do s1 = 1,constitutive_nonlocal_totalNslip(i) + do s2 = 1,constitutive_nonlocal_totalNslip(i) + + constitutive_nonlocal_forestProjectionEdge(s1, s2, i) & + = abs(math_mul3x3(lattice_sn(:, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure), & + lattice_st(:, constitutive_nonlocal_slipSystemLattice(s2,i), myStructure))) ! forest projection of edge dislocations is the projection of (b x n) onto the slip normal of the respective splip plane + + + constitutive_nonlocal_forestProjectionScrew(s1, s2, i) & + = abs(math_mul3x3(lattice_sn(:, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure), & + lattice_sd(:, constitutive_nonlocal_slipSystemLattice(s2,i), myStructure))) ! forest projection of screw dislocations is the projection of b onto the slip normal of the respective splip plane + + enddo; enddo + +!*** calculation of interaction matrices + + do s1 = 1,constitutive_nonlocal_totalNslip(i) + do s2 = 1,constitutive_nonlocal_totalNslip(i) + + constitutive_nonlocal_interactionMatrixSlipSlip(s1, s2, i) & + = constitutive_nonlocal_interactionSlipSlip( lattice_interactionSlipSlip(constitutive_nonlocal_slipSystemLattice(s1,i), & + constitutive_nonlocal_slipSystemLattice(s2,i), & + myStructure), & + i ) + + enddo; enddo + +enddo + +endsubroutine + + + +!********************************************************************* +!* initial microstructural state (just the "basic" states) * +!********************************************************************* +pure function constitutive_nonlocal_stateInit(myInstance) + +use prec, only: pReal, & + pInt +use lattice, only: lattice_maxNslipFamily +use IO, only: IO_error +implicit none + +!*** input variables +integer(pInt), intent(in) :: myInstance ! number specifying the current instance of the constitution + +!*** output variables +real(pReal), dimension(constitutive_nonlocal_sizeState(myInstance)) :: & + constitutive_nonlocal_stateInit + +!*** local variables +real(pReal), dimension(constitutive_nonlocal_totalNslip(myInstance)) :: & + rhoEdgePos, & ! positive edge dislocation density + rhoEdgeNeg, & ! negative edge dislocation density + rhoScrewPos, & ! positive screw dislocation density + rhoScrewNeg, & ! negative screw dislocation density + rhoForest, & ! forest dislocation density + tauSlipThreshold ! threshold shear stress for slip +integer(pInt) ns, & ! short notation for total number of active slip systems + f, & ! index of lattice family + s0, & + s1, & + s ! index of slip system + + +constitutive_nonlocal_stateInit = 0.0_pReal +ns = constitutive_nonlocal_totalNslip(myInstance) + + +!*** set the basic state variables + +s1 = 0_pInt +do f = 1,lattice_maxNslipFamily + + s0 = s1 + 1_pInt + s1 = s0 + constitutive_nonlocal_Nslip(f,myInstance) - 1_pInt + + do s = s0,s1 + rhoEdgePos(s) = constitutive_nonlocal_rhoEdgePos0(f, myInstance) + rhoEdgeNeg(s) = constitutive_nonlocal_rhoEdgeNeg0(f, myInstance) + rhoScrewPos(s) = constitutive_nonlocal_rhoScrewPos0(f, myInstance) + rhoScrewNeg(s) = constitutive_nonlocal_rhoScrewNeg0(f, myInstance) + +enddo; enddo + + +!*** calculate the dependent state variables + +! forest dislocation density +forall (s = 1:ns) & + rhoForest(s) = dot_product( (rhoEdgePos + rhoEdgeNeg), constitutive_nonlocal_forestProjectionEdge(1:ns, s, myInstance) ) & + + dot_product( (rhoScrewPos + rhoScrewNeg), constitutive_nonlocal_forestProjectionScrew(1:ns, s, myInstance) ) ! calculation of forest dislocation density as projection of screw and edge dislocations + + +! threshold shear stress for dislocation slip +forall (s = 1:ns) & + tauSlipThreshold(s) = constitutive_nonlocal_Gmod(myInstance) & + * constitutive_nonlocal_burgersBySlipSystem(s, myInstance) & + * sqrt( dot_product( (rhoEdgePos + rhoEdgeNeg + rhoScrewPos + rhoScrewNeg), & + constitutive_nonlocal_interactionMatrixSlipSlip(1:ns, s, myInstance) ) ) + + +!*** put everything together and in right order + +constitutive_nonlocal_stateInit( 1: ns) = rhoEdgePos +constitutive_nonlocal_stateInit( ns+1:2*ns) = rhoEdgeNeg +constitutive_nonlocal_stateInit(2*ns+1:3*ns) = rhoScrewPos +constitutive_nonlocal_stateInit(3*ns+1:4*ns) = rhoScrewNeg +constitutive_nonlocal_stateInit(4*ns+1:5*ns) = rhoForest +constitutive_nonlocal_stateInit(5*ns+1:6*ns) = tauSlipThreshold + +endfunction + + + +!********************************************************************* +!* calculates homogenized elacticity matrix * +!********************************************************************* +pure function constitutive_nonlocal_homogenizedC(state,g,ip,el) + +use prec, only: pReal, & + pInt, & + p_vec +use mesh, only: mesh_NcpElems, & + mesh_maxNips +use material, only: homogenization_maxNgrains, & + material_phase, & + phase_constitutionInstance +implicit none + +!*** input variables +integer(pInt), intent(in) :: g, & ! current grain ID + ip, & ! current integration point + el ! current element +type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state ! microstructural state + +!*** output variables +real(pReal), dimension(6,6) :: constitutive_nonlocal_homogenizedC ! homogenized elasticity matrix + +!*** local variables +integer(pInt) myInstance ! current instance of this constitution + +myInstance = phase_constitutionInstance(material_phase(g,ip,el)) + +constitutive_nonlocal_homogenizedC = constitutive_nonlocal_Cslip_66(:,:,myInstance) + +endfunction + + + +!********************************************************************* +!* calculates quantities characterizing the microstructure * +!********************************************************************* +subroutine constitutive_nonlocal_microstructure(Temperature, Fp, state, g, ip, el) + +use prec, only: pReal, & + pInt, & + p_vec +use math, only: math_Plain3333to99, & + math_Mandel33to6, & + math_Mandel6to33, & + math_mul33x33, & + math_mul3x3, & + math_mul33x3, & + math_transpose3x3, & + pi +use mesh, only: mesh_NcpElems, & + mesh_maxNips, & + mesh_element, & + FE_NipNeighbors, & + mesh_ipNeighborhood, & + mesh_ipVolume, & + mesh_ipCenterOfGravity +use material, only: homogenization_maxNgrains, & + material_phase, & + phase_constitutionInstance +use lattice, only: lattice_Sslip, & + lattice_Sslip_v, & + lattice_maxNslipFamily, & + lattice_NslipSystem, & + lattice_maxNslip, & + lattice_sd, & + lattice_sn, & + lattice_st + +implicit none + +!*** input variables +integer(pInt), intent(in) :: g, & ! current grain ID + ip, & ! current integration point + el ! current element +real(pReal), intent(in) :: Temperature ! temperature +real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + Fp ! plastic deformation gradient + +!*** input/output variables +type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(inout) :: & + state ! microstructural state + +!*** output variables + +!*** 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 + neighboring_el, & ! element number of my neighbor + neighboring_ip, & ! integration point of my neighbor + n, & ! index of my current neighbor + s, & ! index of my current slip system + sLattice ! index of my current slip system as specified by lattice +real(pReal) gb, & ! short notation for G*b/2/pi + x, & ! coordinate in direction of lvec + y, & ! coordinate in direction of bvec + z ! coordinate in direction of nvec +real(pReal), dimension(3) :: connectingVector ! vector connecting the centers of gravity of me and my neigbor +real(pReal), dimension(6) :: backStress_v ! backstress resulting from the neighboring excess dislocation densities as 2nd Piola-Kirchhoff stress in Mandel notation +real(pReal), dimension(3,3) :: transform, & ! orthogonal transformation matrix from slip coordinate system with e1=bxn, e2=b, e3=n to lattice coordinate system + sigma ! backstress resulting from the excess dislocation density of a single slip system and a single neighbor calculated in the coordinate system of the slip system +real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: & + rhoEdgePos, & ! positive edge dislocation density + rhoEdgeNeg, & ! negative edge dislocation density + rhoScrewPos, & ! positive screw dislocation density + rhoScrewNeg, & ! negative screw dislocation density + rhoForest, & ! forest dislocation density + tauSlipThreshold, & ! threshold shear stress + neighboring_rhoEdgePos, & ! positive edge dislocation density of my neighbor + neighboring_rhoEdgeNeg, & ! negative edge dislocation density of my neighbor + neighboring_rhoScrewPos, & ! positive screw dislocation density of my neighbor + neighboring_rhoScrewNeg, & ! negative screw dislocation density of my neighbor + neighboring_rhoEdgeExcess, &! edge excess dislocation density of my neighbor + neighboring_rhoScrewExcess,&! screw excess dislocation density of my neighbor + neighboring_Nedge, & ! total number of edge excess dislocations in my neighbor + neighboring_Nscrew + + +myInstance = phase_constitutionInstance(material_phase(g,ip,el)) +myStructure = constitutive_nonlocal_structure(myInstance) +ns = constitutive_nonlocal_totalNslip(myInstance) + + +!********************************************************************** +!*** get basic states + +rhoEdgePos = state(g,ip,el)%p( 1: ns) +rhoEdgeNeg = state(g,ip,el)%p( ns+1:2*ns) +rhoScrewPos = state(g,ip,el)%p(2*ns+1:3*ns) +rhoScrewNeg = state(g,ip,el)%p(3*ns+1:4*ns) + + +!********************************************************************** +!*** calculate dependent states + +!*** calculate the forest dislocation density + +forall (s = 1:ns) & + rhoForest(s) = dot_product( (rhoEdgePos + rhoEdgeNeg), constitutive_nonlocal_forestProjectionEdge(1:ns, s, myInstance) ) & + + dot_product( (rhoScrewPos + rhoScrewNeg), constitutive_nonlocal_forestProjectionScrew(1:ns, s, myInstance) ) ! calculation of forest dislocation density as projection of screw and edge dislocations + + +!*** calculate the threshold shear stress for dislocation slip + +forall (s = 1:ns) & + tauSlipThreshold(s) = constitutive_nonlocal_Gmod(myInstance) & + * constitutive_nonlocal_burgersBySlipSystem(s, myInstance) & + * sqrt( dot_product( (rhoEdgePos + rhoEdgeNeg + rhoScrewPos + rhoScrewNeg), & + constitutive_nonlocal_interactionMatrixSlipSlip(1:ns, s, myInstance) ) ) + + +!*** calculate the backstress of the neighboring excess dislocation densities + +backStress_v = 0.0_pReal + + +! loop through my neighbors (if it exists!) + +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) + + if ( neighboring_el == 0 .or. neighboring_ip == 0 ) cycle + + + ! calculate the connecting vector between me and my neighbor and his excess dislocation density + + connectingVector = math_mul33x3( Fp(:,:,g,ip,el), & + (mesh_ipCenterOfGravity(:,ip,el) - mesh_ipCenterOfGravity(:,neighboring_ip,neighboring_el)) ) + + neighboring_rhoEdgePos = state(1, neighboring_ip, neighboring_el)%p( 1: ns) + neighboring_rhoEdgeNeg = state(1, neighboring_ip, neighboring_el)%p( ns+1:2*ns) + neighboring_rhoScrewPos = state(1, neighboring_ip, neighboring_el)%p(2*ns+1:3*ns) + neighboring_rhoScrewNeg = state(1, neighboring_ip, neighboring_el)%p(3*ns+1:4*ns) + + neighboring_rhoEdgeExcess = neighboring_rhoEdgePos - neighboring_rhoEdgeNeg + neighboring_rhoScrewExcess = neighboring_rhoScrewPos - neighboring_rhoScrewNeg + + neighboring_Nedge = neighboring_rhoEdgeExcess * mesh_ipVolume(neighboring_ip, neighboring_el) ** (2.0_pReal/3.0_pReal) + neighboring_Nscrew = neighboring_rhoScrewExcess * mesh_ipVolume(neighboring_ip, neighboring_el) ** (2.0_pReal/3.0_pReal) + + ! loop over slip systems and get their slip directions, slip normals, and sd x sn + + do s = 1,ns + + transform = reshape( (/lattice_st(:, constitutive_nonlocal_slipSystemLattice(s,myInstance), myStructure), & + lattice_sd(:, constitutive_nonlocal_slipSystemLattice(s,myInstance), myStructure), & + lattice_sn(:, constitutive_nonlocal_slipSystemLattice(s,myInstance), myStructure)/), (/3,3/) ) + + + ! coordinate transformation of p from the lattice coordinate system to the slip coordinate system + x = math_mul3x3(connectingVector, transform(:,1)) + y = math_mul3x3(connectingVector, transform(:,2)) + z = math_mul3x3(connectingVector, transform(:,3)) + + + ! calculate the back stress in the slip coordinate system for this slip system + gb = constitutive_nonlocal_Gmod(myInstance) * constitutive_nonlocal_burgersBySlipSystem(s,myInstance) / (2.0_pReal*pi) + + sigma(2,2) = - gb * neighboring_Nedge(s) / (1.0_pReal-constitutive_nonlocal_nu(myInstance)) & + * z * (3.0_pReal*y**2.0_pReal + z**2.0_pReal) / (y**2.0_pReal + z**2.0_pReal)**2.0_pReal + + sigma(3,3) = gb * neighboring_Nedge(s) / (1.0_pReal-constitutive_nonlocal_nu(myInstance)) & + * z * (y**2.0_pReal - z**2.0_pReal) / (y**2.0_pReal + z**2.0_pReal)**2.0_pReal + + sigma(1,1) = constitutive_nonlocal_nu(myInstance) * (sigma(2,2) + sigma(3,3)) + + sigma(1,2) = gb * neighboring_Nscrew(s) * z / (x**2.0_pReal + z**2.0_pReal) + + sigma(2,3) = gb * ( neighboring_Nedge(s) / (1.0_pReal-constitutive_nonlocal_nu(myInstance)) & + * (y**2.0_pReal - z**2.0_pReal) / (y**2.0_pReal + z**2.0_pReal) & + - neighboring_Nscrew(s) * x / (x**2.0_pReal + z**2.0_pReal) ) + + sigma(2,1) = sigma(1,2) + sigma(3,2) = sigma(2,3) + sigma(1,3) = 0.0_pReal + sigma(3,1) = 0.0_pReal + + ! coordinate transformation from the slip coordinate system to the lattice coordinate system + backStress_v = backStress_v + math_Mandel33to6( math_mul33x33(math_transpose3x3(transform), math_mul33x33(sigma, transform) ) ) + + enddo +enddo + + +!********************************************************************** +!*** set dependent states + +state(g,ip,el)%p(4*ns+1:5*ns) = rhoForest +state(g,ip,el)%p(5*ns+1:6*ns) = tauSlipThreshold +state(g,ip,el)%p(6*ns+1:6*ns+6) = backstress_v + +endsubroutine + + + +!********************************************************************* +!* calculates plastic velocity gradient and its tangent * +!********************************************************************* +subroutine constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar99, Tstar_v, Temperature, state, g, ip, el) + +use prec, only: pReal, & + pInt, & + p_vec +use math, only: math_Plain3333to99, & + math_mul6x6, & + math_Mandel6to33 +use debug, only: debugger +use mesh, only: mesh_NcpElems, & + mesh_maxNips +use material, only: homogenization_maxNgrains, & + material_phase, & + phase_constitutionInstance +use lattice, only: lattice_Sslip, & + lattice_Sslip_v + +implicit none + +!*** input variables +integer(pInt), intent(in) :: g, & ! current grain number + ip, & ! current integration point + el ! current element number +real(pReal), intent(in) :: Temperature ! temperature +type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + state ! microstructural state +real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola-Kirchhoff stress in Mandel notation + +!*** output variables +real(pReal), dimension(3,3), intent(out) :: Lp ! plastic velocity gradient +real(pReal), dimension(9,9), intent(out) :: dLp_dTstar99 ! derivative of Lp with respect to Tstar (9x9 matrix) + +!*** 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 + i, & + j, & + k, & + l, & + t, & ! dislocation type + s, & ! index of my current slip system + sLattice ! index of my current slip system as specified by lattice +real(pReal), dimension(6) :: backStress_v ! backstress resulting from the neighboring excess dislocation densities as 2nd Piola-Kirchhoff stress +real(pReal), dimension(3,3,3,3) :: dLp_dTstar3333 ! derivative of Lp with respect to Tstar (3x3x3x3 matrix) +real(pReal), dimension(4,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: & + rho ! dislocation densities +real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: & + rhoForest, & ! forest dislocation density + tauSlipThreshold, & ! threshold shear stress + tauSlip, & ! resolved shear stress + gdotSlip, & ! shear rate + dgdot_dtauSlip, & ! derivative of the shear rate with respect to the shear stress + v ! dislocation velocity + + +!*** initialize local variables + +v = 0.0_pReal +tauSlip = 0.0_pReal +gdotSlip = 0.0_pReal +Lp = 0.0_pReal +dLp_dTstar3333 = 0.0_pReal + +myInstance = phase_constitutionInstance(material_phase(g,ip,el)) +myStructure = constitutive_nonlocal_structure(myInstance) +ns = constitutive_nonlocal_totalNslip(myInstance) + +!*** shortcut to state variables + +forall (t = 1:4) rho(t,:) = state(g,ip,el)%p((t-1)*ns+1:t*ns) +rhoForest = state(g,ip,el)%p(4*ns+1:5*ns) +tauSlipThreshold = state(g,ip,el)%p(5*ns+1:6*ns) +backStress_v = state(g,ip,el)%p(6*ns+1:6*ns+6) +! if (debugger) write(6,'(a20,3(i3,x),/,3(3(f12.3,x)/))') 'backstress / MPa at ', g,ip,el, math_Mandel6to33(backStress_v/1e6) +! if (debugger) write(6,'(a15,3(i3,x),/,3(3(f12.3,x)/))') 'Tstar / MPa at ',g,ip,el, math_Mandel6to33(Tstar_v/1e6) + +!*** loop over slip systems + +do s = 1,constitutive_nonlocal_totalNslip(myInstance) + + sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance) + + + !*** Calculation of Lp + + tauSlip(s) = math_mul6x6( Tstar_v + backStress_v, lattice_Sslip_v(:,sLattice,myStructure) ) + + if (rhoForest(s) > 0.0_pReal) & + v(s) = constitutive_nonlocal_v0BySlipSystem(s,myInstance) & + * exp( - ( tauSlipThreshold(s) - abs(tauSlip(s)) ) * constitutive_nonlocal_burgersBySlipSystem(s,myInstance)**2.0_pReal & + / ( kB * Temperature * sqrt(rhoForest(s)) ) ) & + * sign(1.0_pReal,tauSlip(s)) + + gdotSlip(s) = sum(rho(:,s)) * constitutive_nonlocal_burgersBySlipSystem(s,myInstance) * v(s) + + Lp = Lp + gdotSlip(s) * lattice_Sslip(:,:,sLattice,myStructure) + + + !*** Calculation of the tangent of Lp + + dgdot_dtauSlip(s) = gdotSlip(s) * constitutive_nonlocal_burgersBySlipSystem(s,myInstance)**2.0_pReal & + / ( kB * Temperature * sqrt(rhoForest(s)) ) + + forall (i=1:3,j=1:3,k=1:3,l=1:3) & + dLp_dTstar3333(i,j,k,l) = dLp_dTstar3333(i,j,k,l) + dgdot_dtauSlip(s) * lattice_Sslip(i,j, sLattice,myStructure) & + * lattice_Sslip(k,l, sLattice,myStructure) +enddo +!if (debugger) write(6,'(a26,3(i3,x),/,12(f10.3,x),/)') 'tauSlipThreshold / MPa at ',g,ip,el, tauSlipThreshold/1e6 +!if (debugger) write(6,'(a15,3(i3,x),/,12(f10.3,x),/)') 'tauSlip / MPa at ',g,ip,el, tauSlip/1e6 +!if (debugger) write(6,'(a15,3(i3,x),/,12(e10.3,x),/)') 'gdotSlip at ',g,ip,el, gdotSlip + +dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333) + +endsubroutine + + + +!********************************************************************* +!* rate of change of microstructure * +!********************************************************************* +subroutine constitutive_nonlocal_dotState(dotState, Tstar_v, Fp, invFp, Temperature, state, g, ip, el) + +use prec, only: pReal, & + pInt, & + p_vec +use debug, only: debugger +use math, only: math_norm3, & + math_mul6x6, & + math_mul3x3, & + math_mul33x3, & + math_transpose3x3 +use mesh, only: mesh_NcpElems, & + mesh_maxNips, & + mesh_element, & + FE_NipNeighbors, & + mesh_ipNeighborhood, & + mesh_ipVolume, & + mesh_ipCenterOfGravity, & + mesh_ipArea, & + mesh_ipAreaNormal +use material, only: homogenization_maxNgrains, & + material_phase, & + phase_constitutionInstance +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) :: Temperature ! temperature +real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola-Kirchhoff stress in Mandel notation +real(pReal), dimension(3,3), intent(in) :: Fp, & ! plastic deformation gradient + invFp ! inverse of plastic deformation gradient +type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + state ! microstructural state + +!*** input/output variables +type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(inout) :: & + dotState ! evolution of state variables / microstructure + +!*** output variables + +!*** 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 + neighboring_el, & ! element number of my neighbor + neighboring_ip, & ! integration point of my neighbor + n, & ! index of my current neighbor + t, & ! type of dislocation + s, & ! index of my current slip system + sLattice ! index of my current slip system as specified by lattice +real(pReal), dimension(4,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: & + rho, & ! dislocation densities + gdot ! shear rates +real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: & + rhoForest, & ! forest dislocation density + tauSlipThreshold, & ! threshold shear stress + tauSlip, & ! resolved shear stress + v ! dislocation velocity +real(pReal), dimension(3,4,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: & + m ! direction of dislocation motion +real(pReal), dimension(6) :: backStress_v ! backstress resulting from the neighboring excess dislocation densities as 2nd Piola-Kirchhoff stress +real(pReal), dimension(4) :: leavingRho ! dislocation densities leaving the current interface +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 + + + +myInstance = phase_constitutionInstance(material_phase(g,ip,el)) +myStructure = constitutive_nonlocal_structure(myInstance) +ns = constitutive_nonlocal_totalNslip(myInstance) + + +!*** shortcut to state variables + +forall (t = 1:4) rho(t,:) = state(g,ip,el)%p((t-1)*ns+1:t*ns) +rhoForest = state(g,ip,el)%p(4*ns+1:5*ns) +tauSlipThreshold = state(g,ip,el)%p(5*ns+1:6*ns) +backStress_v = state(g,ip,el)%p(6*ns+1:6*ns+6) + + +do s = 1,ns + +!*** Calculate gdot for each type of dislocation + + sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance) + tauSlip(s) = math_mul6x6( Tstar_v + backStress_v, lattice_Sslip_v(:,sLattice,myStructure) ) + + v(s) = constitutive_nonlocal_v0BySlipSystem(s,myInstance) & + * exp( - ( tauSlipThreshold(s) - abs(tauSlip(s)) ) & + * constitutive_nonlocal_burgersBySlipSystem(s,myInstance)**2.0_pReal & + / ( kB * Temperature * sqrt(rhoForest(s)) ) ) & + * sign(1.0_pReal,tauSlip(s)) + + forall (t = 1:4) gdot(t,s) = rho(t,s) * constitutive_nonlocal_burgersBySlipSystem(s,myInstance) * v(s) + + +!*** Direction of dislocation motion + + m(:,1,s) = lattice_sd(:, sLattice, myStructure) + m(:,2,s) = -lattice_sd(:, sLattice, myStructure) + m(:,3,s) = lattice_st(:, sLattice, myStructure) + m(:,4,s) = -lattice_st(:, sLattice, myStructure) + +enddo + + +!*** loop through my neighbors if existent and calculate the area and the surface normal of the interface + +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) + + surfaceNormal = math_mul33x3(math_transpose3x3(invFp), mesh_ipAreaNormal(:,n,ip,el)) + norm_surfaceNormal = math_norm3(surfaceNormal) + surfaceNormal = surfaceNormal / norm_surfaceNormal + area = mesh_ipArea(n,ip,el) / norm_surfaceNormal + + +!*** loop through my interfaces +!*** subtract dislocation densities that leave through an interface from my dotState and add them to the neighboring dotState + + do s = 1,ns + + leavingRho = 0.0_pReal + do t = 1,4 + if ( sign(1.0_pReal,math_mul3x3(m(:,t,s),surfaceNormal)) == sign(1.0_pReal,gdot(t,s)) ) then + + leavingRho(t) = gdot(t,s) / constitutive_nonlocal_burgersBySlipSystem(s,myInstance) & + * math_mul3x3(m(:,t,s),surfaceNormal) * area + + dotState(1,ip,el)%p((t-1)*ns+s) = dotState(1,ip,el)%p((t-1)*ns+s) - leavingRho(t) + + 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) + leavingRho(t) + endif + endif + enddo + + enddo + +enddo + +endsubroutine + + + +!********************************************************************* +!* rate of change of temperature * +!********************************************************************* +pure function constitutive_nonlocal_dotTemperature(Tstar_v,Temperature,state,g,ip,el) + +use prec, only: pReal, & + pInt, & + p_vec +use mesh, only: mesh_NcpElems, & + mesh_maxNips +use material, only: homogenization_maxNgrains +implicit none + +!* input variables +integer(pInt), intent(in) :: g, & ! current grain ID + ip, & ! current integration point + el ! current element +real(pReal), intent(in) :: Temperature ! temperature +real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola-Kirchhoff stress in Mandel notation +type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + state ! microstructural state + +!* output variables +real(pReal) constitutive_nonlocal_dotTemperature ! evolution of Temperature + +!* local variables + +constitutive_nonlocal_dotTemperature = 0.0_pReal + +endfunction + + + +!********************************************************************* +!* return array of constitutive results * +!* INPUT: * +!* - Temperature : temperature * +!* - Tstar_v : 2nd Piola Kirchhoff stress tensor (Mandel) * +!* - dt : current time increment * +!* - ipc : component-ID at current integration point * +!* - ip : current integration point * +!* - el : current element * +!********************************************************************* +pure function constitutive_nonlocal_postResults(Tstar_v,Temperature,dt,state,ipc,ip,el) + +use prec, only: pReal,pInt,p_vec +use mesh, only: mesh_NcpElems,mesh_maxNips +use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance,phase_Noutput +use lattice, only: lattice_Sslip_v,lattice_Stwin_v,lattice_maxNslipFamily,lattice_maxNtwinFamily, & + lattice_NslipSystem,lattice_NtwinSystem,lattice_shearTwin +implicit none + +!* Definition of variables +integer(pInt), intent(in) :: ipc,ip,el +real(pReal), intent(in) :: dt,Temperature +real(pReal), dimension(6), intent(in) :: Tstar_v +type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state +integer(pInt) matID,structID,ns,nt,f,o,i,c,j,index_myFamily +real(pReal) sumf,tau +real(pReal), dimension(constitutive_nonlocal_sizePostResults(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: & + constitutive_nonlocal_postResults + +!* Shortened notation +matID = phase_constitutionInstance(material_phase(ipc,ip,el)) +structID = constitutive_nonlocal_structure(matID) +ns = constitutive_nonlocal_totalNslip(matID) + + +!* Required output +c = 0_pInt +constitutive_nonlocal_postResults = 0.0_pReal + +do o = 1,phase_Noutput(material_phase(ipc,ip,el)) + select case(constitutive_nonlocal_output(o,matID)) + + case ('dislocationdensity') + constitutive_nonlocal_postResults(c+1:c+ns) = state(ipc,ip,el)%p(1:ns) + c = c + ns + + end select +enddo + +endfunction + +END MODULE \ No newline at end of file diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 0ff6162c3..f69226f3e 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -34,6 +34,7 @@ real(pReal), dimension (:,:,:,:), allocatable :: crystallite_Tstar_v, & crystallite_subTstar0_v ! 2nd Piola-Kirchhoff stress vector at start of crystallite inc real(pReal), dimension (:,:,:,:,:), allocatable :: crystallite_Fe, & ! current "elastic" def grad (end of converged time step) crystallite_Fp, & ! current plastic def grad (end of converged time step) + crystallite_invFp, & ! inverse of current plastic def grad (end of converged time step) crystallite_Fp0, & ! plastic def grad at start of FE inc crystallite_partionedFp0,& ! plastic def grad at start of homog inc crystallite_subFp0,& ! plastic def grad at start of crystallite inc @@ -53,7 +54,10 @@ real(pReal), dimension (:,:,:,:,:,:,:), allocatable :: crystallite_dPdF, & logical, dimension (:,:,:), allocatable :: crystallite_localConstitution, & ! indicates this grain to have purely local constitutive law crystallite_requested, & ! flag to request crystallite calculation crystallite_onTrack, & ! flag to indicate ongoing calculation - crystallite_converged ! convergence flag + crystallite_converged, & ! convergence flag + crystallite_stateConverged, & ! flag indicating convergence of state + crystallite_temperatureConverged, & ! flag indicating convergence of temperature + crystallite_nonfinished ! requested and ontrack but not converged CONTAINS @@ -99,38 +103,42 @@ subroutine crystallite_init(Temperature) iMax = mesh_maxNips eMax = mesh_NcpElems - allocate(crystallite_Temperature(gMax,iMax,eMax)); crystallite_Temperature = Temperature - allocate(crystallite_P(3,3,gMax,iMax,eMax)); crystallite_P = 0.0_pReal - allocate(crystallite_Fe(3,3,gMax,iMax,eMax)); crystallite_Fe = 0.0_pReal - allocate(crystallite_Fp(3,3,gMax,iMax,eMax)); crystallite_Fp = 0.0_pReal - allocate(crystallite_Lp(3,3,gMax,iMax,eMax)); crystallite_Lp = 0.0_pReal - allocate(crystallite_Tstar_v(6,gMax,iMax,eMax)); crystallite_Tstar_v = 0.0_pReal - allocate(crystallite_F0(3,3,gMax,iMax,eMax)); crystallite_F0 = 0.0_pReal - allocate(crystallite_Fp0(3,3,gMax,iMax,eMax)); crystallite_Fp0 = 0.0_pReal - allocate(crystallite_Lp0(3,3,gMax,iMax,eMax)); crystallite_Lp0 = 0.0_pReal - allocate(crystallite_Tstar0_v(6,gMax,iMax,eMax)); crystallite_Tstar0_v = 0.0_pReal - allocate(crystallite_partionedTemperature0(gMax,iMax,eMax)); crystallite_partionedTemperature0 = 0.0_pReal - allocate(crystallite_partionedF(3,3,gMax,iMax,eMax)); crystallite_partionedF = 0.0_pReal - allocate(crystallite_partionedF0(3,3,gMax,iMax,eMax)); crystallite_partionedF0 = 0.0_pReal - allocate(crystallite_partionedFp0(3,3,gMax,iMax,eMax)); crystallite_partionedFp0 = 0.0_pReal - allocate(crystallite_partionedLp0(3,3,gMax,iMax,eMax)); crystallite_partionedLp0 = 0.0_pReal - allocate(crystallite_partionedTstar0_v(6,gMax,iMax,eMax)); crystallite_partionedTstar0_v = 0.0_pReal - allocate(crystallite_subTemperature0(gMax,iMax,eMax)); crystallite_subTemperature0 = 0.0_pReal - allocate(crystallite_subF(3,3,gMax,iMax,eMax)); crystallite_subF = 0.0_pReal - allocate(crystallite_subF0(3,3,gMax,iMax,eMax)); crystallite_subF0 = 0.0_pReal - allocate(crystallite_subFp0(3,3,gMax,iMax,eMax)); crystallite_subFp0 = 0.0_pReal - allocate(crystallite_subLp0(3,3,gMax,iMax,eMax)); crystallite_subLp0 = 0.0_pReal - allocate(crystallite_subTstar0_v(6,gMax,iMax,eMax)); crystallite_subTstar0_v = 0.0_pReal - allocate(crystallite_dPdF(3,3,3,3,gMax,iMax,eMax)); crystallite_dPdF = 0.0_pReal - allocate(crystallite_fallbackdPdF(3,3,3,3,gMax,iMax,eMax)); crystallite_fallbackdPdF = 0.0_pReal - allocate(crystallite_dt(gMax,iMax,eMax)); crystallite_dt = 0.0_pReal - allocate(crystallite_subdt(gMax,iMax,eMax)); crystallite_subdt = 0.0_pReal - allocate(crystallite_subFrac(gMax,iMax,eMax)); crystallite_subFrac = 0.0_pReal - allocate(crystallite_subStep(gMax,iMax,eMax)); crystallite_subStep = 0.0_pReal - allocate(crystallite_localConstitution(gMax,iMax,eMax)); crystallite_localConstitution = .true. - allocate(crystallite_requested(gMax,iMax,eMax)); crystallite_requested = .false. - allocate(crystallite_onTrack(gMax,iMax,eMax)); crystallite_onTrack = .false. - allocate(crystallite_converged(gMax,iMax,eMax)); crystallite_converged = .true. + allocate(crystallite_Temperature(gMax,iMax,eMax)); crystallite_Temperature = Temperature + allocate(crystallite_P(3,3,gMax,iMax,eMax)); crystallite_P = 0.0_pReal + allocate(crystallite_Fe(3,3,gMax,iMax,eMax)); crystallite_Fe = 0.0_pReal + allocate(crystallite_Fp(3,3,gMax,iMax,eMax)); crystallite_Fp = 0.0_pReal + allocate(crystallite_invFp(3,3,gMax,iMax,eMax)); crystallite_invFp = 0.0_pReal + allocate(crystallite_Lp(3,3,gMax,iMax,eMax)); crystallite_Lp = 0.0_pReal + allocate(crystallite_Tstar_v(6,gMax,iMax,eMax)); crystallite_Tstar_v = 0.0_pReal + allocate(crystallite_F0(3,3,gMax,iMax,eMax)); crystallite_F0 = 0.0_pReal + allocate(crystallite_Fp0(3,3,gMax,iMax,eMax)); crystallite_Fp0 = 0.0_pReal + allocate(crystallite_Lp0(3,3,gMax,iMax,eMax)); crystallite_Lp0 = 0.0_pReal + allocate(crystallite_Tstar0_v(6,gMax,iMax,eMax)); crystallite_Tstar0_v = 0.0_pReal + allocate(crystallite_partionedTemperature0(gMax,iMax,eMax)); crystallite_partionedTemperature0 = 0.0_pReal + allocate(crystallite_partionedF(3,3,gMax,iMax,eMax)); crystallite_partionedF = 0.0_pReal + allocate(crystallite_partionedF0(3,3,gMax,iMax,eMax)); crystallite_partionedF0 = 0.0_pReal + allocate(crystallite_partionedFp0(3,3,gMax,iMax,eMax)); crystallite_partionedFp0 = 0.0_pReal + allocate(crystallite_partionedLp0(3,3,gMax,iMax,eMax)); crystallite_partionedLp0 = 0.0_pReal + allocate(crystallite_partionedTstar0_v(6,gMax,iMax,eMax)); crystallite_partionedTstar0_v = 0.0_pReal + allocate(crystallite_subTemperature0(gMax,iMax,eMax)); crystallite_subTemperature0 = 0.0_pReal + allocate(crystallite_subF(3,3,gMax,iMax,eMax)); crystallite_subF = 0.0_pReal + allocate(crystallite_subF0(3,3,gMax,iMax,eMax)); crystallite_subF0 = 0.0_pReal + allocate(crystallite_subFp0(3,3,gMax,iMax,eMax)); crystallite_subFp0 = 0.0_pReal + allocate(crystallite_subLp0(3,3,gMax,iMax,eMax)); crystallite_subLp0 = 0.0_pReal + allocate(crystallite_subTstar0_v(6,gMax,iMax,eMax)); crystallite_subTstar0_v = 0.0_pReal + allocate(crystallite_dPdF(3,3,3,3,gMax,iMax,eMax)); crystallite_dPdF = 0.0_pReal + allocate(crystallite_fallbackdPdF(3,3,3,3,gMax,iMax,eMax)); crystallite_fallbackdPdF = 0.0_pReal + allocate(crystallite_dt(gMax,iMax,eMax)); crystallite_dt = 0.0_pReal + allocate(crystallite_subdt(gMax,iMax,eMax)); crystallite_subdt = 0.0_pReal + allocate(crystallite_subFrac(gMax,iMax,eMax)); crystallite_subFrac = 0.0_pReal + allocate(crystallite_subStep(gMax,iMax,eMax)); crystallite_subStep = 0.0_pReal + allocate(crystallite_localConstitution(gMax,iMax,eMax)); crystallite_localConstitution = .true. + allocate(crystallite_requested(gMax,iMax,eMax)); crystallite_requested = .false. + allocate(crystallite_onTrack(gMax,iMax,eMax)); crystallite_onTrack = .true. + allocate(crystallite_converged(gMax,iMax,eMax)); crystallite_converged = .true. + allocate(crystallite_stateConverged(gMax,iMax,eMax)); crystallite_stateConverged = .false. + allocate(crystallite_temperatureConverged(gMax,iMax,eMax)); crystallite_temperatureConverged = .false. + allocate(crystallite_nonfinished(gMax,iMax,eMax)); crystallite_nonfinished = .true. !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over all cp elements @@ -149,7 +157,7 @@ subroutine crystallite_init(Temperature) enddo enddo !$OMPEND PARALLEL DO - + call crystallite_stressAndItsTangent(.true.) ! request elastic answers crystallite_fallbackdPdF = crystallite_dPdF ! use initial elastic stiffness as fallback @@ -158,39 +166,43 @@ subroutine crystallite_init(Temperature) write(6,*) write(6,*) '<<<+- crystallite init -+>>>' write(6,*) - write(6,'(a32,x,7(i5,x))') 'crystallite_Nresults: ', crystallite_Nresults - write(6,'(a32,x,7(i5,x))') 'crystallite_Temperature ', shape(crystallite_Temperature) - write(6,'(a32,x,7(i5,x))') 'crystallite_Fe: ', shape(crystallite_Fe) - write(6,'(a32,x,7(i5,x))') 'crystallite_Fp: ', shape(crystallite_Fp) - write(6,'(a32,x,7(i5,x))') 'crystallite_Lp: ', shape(crystallite_Lp) - write(6,'(a32,x,7(i5,x))') 'crystallite_F0: ', shape(crystallite_F0) - write(6,'(a32,x,7(i5,x))') 'crystallite_Fp0: ', shape(crystallite_Fp0) - write(6,'(a32,x,7(i5,x))') 'crystallite_Lp0: ', shape(crystallite_Lp0) - write(6,'(a32,x,7(i5,x))') 'crystallite_partionedF: ', shape(crystallite_partionedF) - write(6,'(a32,x,7(i5,x))') 'crystallite_partionedTemp0: ', shape(crystallite_partionedTemperature0) - write(6,'(a32,x,7(i5,x))') 'crystallite_partionedF0: ', shape(crystallite_partionedF0) - write(6,'(a32,x,7(i5,x))') 'crystallite_partionedFp0: ', shape(crystallite_partionedFp0) - write(6,'(a32,x,7(i5,x))') 'crystallite_partionedLp0: ', shape(crystallite_partionedLp0) - write(6,'(a32,x,7(i5,x))') 'crystallite_subF: ', shape(crystallite_subF) - write(6,'(a32,x,7(i5,x))') 'crystallite_subTemperature0: ', shape(crystallite_subTemperature0) - write(6,'(a32,x,7(i5,x))') 'crystallite_subF0: ', shape(crystallite_subF0) - write(6,'(a32,x,7(i5,x))') 'crystallite_subFp0: ', shape(crystallite_subFp0) - write(6,'(a32,x,7(i5,x))') 'crystallite_subLp0: ', shape(crystallite_subLp0) - write(6,'(a32,x,7(i5,x))') 'crystallite_P: ', shape(crystallite_P) - write(6,'(a32,x,7(i5,x))') 'crystallite_Tstar_v: ', shape(crystallite_Tstar_v) - write(6,'(a32,x,7(i5,x))') 'crystallite_Tstar0_v: ', shape(crystallite_Tstar0_v) - write(6,'(a32,x,7(i5,x))') 'crystallite_partionedTstar0_v: ', shape(crystallite_partionedTstar0_v) - write(6,'(a32,x,7(i5,x))') 'crystallite_subTstar0_v: ', shape(crystallite_subTstar0_v) - write(6,'(a32,x,7(i5,x))') 'crystallite_dPdF: ', shape(crystallite_dPdF) - write(6,'(a32,x,7(i5,x))') 'crystallite_fallbackdPdF: ', shape(crystallite_fallbackdPdF) - write(6,'(a32,x,7(i5,x))') 'crystallite_dt: ', shape(crystallite_dt) - write(6,'(a32,x,7(i5,x))') 'crystallite_subdt: ', shape(crystallite_subdt) - write(6,'(a32,x,7(i5,x))') 'crystallite_subFrac: ', shape(crystallite_subFrac) - write(6,'(a32,x,7(i5,x))') 'crystallite_subStep: ', shape(crystallite_subStep) - write(6,'(a32,x,7(i5,x))') 'crystallite_localConstitution: ', shape(crystallite_localConstitution) - write(6,'(a32,x,7(i5,x))') 'crystallite_requested: ', shape(crystallite_requested) - write(6,'(a32,x,7(i5,x))') 'crystallite_onTrack: ', shape(crystallite_onTrack) - write(6,'(a32,x,7(i5,x))') 'crystallite_converged: ', shape(crystallite_converged) + write(6,'(a35,x,7(i5,x))') 'crystallite_Nresults: ', crystallite_Nresults + write(6,*) + write(6,'(a35,x,7(i5,x))') 'crystallite_Temperature: ', shape(crystallite_Temperature) + write(6,'(a35,x,7(i5,x))') 'crystallite_Fe: ', shape(crystallite_Fe) + write(6,'(a35,x,7(i5,x))') 'crystallite_Fp: ', shape(crystallite_Fp) + write(6,'(a35,x,7(i5,x))') 'crystallite_Lp: ', shape(crystallite_Lp) + write(6,'(a35,x,7(i5,x))') 'crystallite_F0: ', shape(crystallite_F0) + write(6,'(a35,x,7(i5,x))') 'crystallite_Fp0: ', shape(crystallite_Fp0) + write(6,'(a35,x,7(i5,x))') 'crystallite_Lp0: ', shape(crystallite_Lp0) + write(6,'(a35,x,7(i5,x))') 'crystallite_partionedF: ', shape(crystallite_partionedF) + write(6,'(a35,x,7(i5,x))') 'crystallite_partionedTemp0: ', shape(crystallite_partionedTemperature0) + write(6,'(a35,x,7(i5,x))') 'crystallite_partionedF0: ', shape(crystallite_partionedF0) + write(6,'(a35,x,7(i5,x))') 'crystallite_partionedFp0: ', shape(crystallite_partionedFp0) + write(6,'(a35,x,7(i5,x))') 'crystallite_partionedLp0: ', shape(crystallite_partionedLp0) + write(6,'(a35,x,7(i5,x))') 'crystallite_subF: ', shape(crystallite_subF) + write(6,'(a35,x,7(i5,x))') 'crystallite_subTemperature0: ', shape(crystallite_subTemperature0) + write(6,'(a35,x,7(i5,x))') 'crystallite_subF0: ', shape(crystallite_subF0) + write(6,'(a35,x,7(i5,x))') 'crystallite_subFp0: ', shape(crystallite_subFp0) + write(6,'(a35,x,7(i5,x))') 'crystallite_subLp0: ', shape(crystallite_subLp0) + write(6,'(a35,x,7(i5,x))') 'crystallite_P: ', shape(crystallite_P) + write(6,'(a35,x,7(i5,x))') 'crystallite_Tstar_v: ', shape(crystallite_Tstar_v) + write(6,'(a35,x,7(i5,x))') 'crystallite_Tstar0_v: ', shape(crystallite_Tstar0_v) + write(6,'(a35,x,7(i5,x))') 'crystallite_partionedTstar0_v: ', shape(crystallite_partionedTstar0_v) + write(6,'(a35,x,7(i5,x))') 'crystallite_subTstar0_v: ', shape(crystallite_subTstar0_v) + write(6,'(a35,x,7(i5,x))') 'crystallite_dPdF: ', shape(crystallite_dPdF) + write(6,'(a35,x,7(i5,x))') 'crystallite_fallbackdPdF: ', shape(crystallite_fallbackdPdF) + write(6,'(a35,x,7(i5,x))') 'crystallite_dt: ', shape(crystallite_dt) + write(6,'(a35,x,7(i5,x))') 'crystallite_subdt: ', shape(crystallite_subdt) + write(6,'(a35,x,7(i5,x))') 'crystallite_subFrac: ', shape(crystallite_subFrac) + write(6,'(a35,x,7(i5,x))') 'crystallite_subStep: ', shape(crystallite_subStep) + write(6,'(a35,x,7(i5,x))') 'crystallite_localConstitution: ', shape(crystallite_localConstitution) + write(6,'(a35,x,7(i5,x))') 'crystallite_requested: ', shape(crystallite_requested) + write(6,'(a35,x,7(i5,x))') 'crystallite_onTrack: ', shape(crystallite_onTrack) + write(6,'(a35,x,7(i5,x))') 'crystallite_converged: ', shape(crystallite_converged) + write(6,'(a35,x,7(i5,x))') 'crystallite_stateConverged: ', shape(crystallite_stateConverged) + write(6,'(a35,x,7(i5,x))') 'crystallite_temperatureConverged: ', shape(crystallite_temperatureConverged) + write(6,'(a35,x,7(i5,x))') 'crystallite_nonfinished: ', shape(crystallite_nonfinished) write(6,*) write(6,*) 'Number of non-local grains: ',count(.not. crystallite_localConstitution) call flush(6) @@ -219,7 +231,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) nCryst use debug, only: debugger, & debug_CrystalliteLoopDistribution, & - debug_StateLoopDistribution, & + debug_CrystalliteStateLoopDistribution, & debug_StiffnessStateLoopDistribution use IO, only: IO_warning use math, only: math_inv3x3, & @@ -235,10 +247,14 @@ subroutine crystallite_stressAndItsTangent(updateJaco) use material, only: homogenization_Ngrains use constitutive, only: constitutive_maxSizeState, & constitutive_sizeState, & + constitutive_sizeDotState, & constitutive_state, & constitutive_subState0, & constitutive_partionedState0, & - constitutive_homogenizedC + constitutive_homogenizedC, & + constitutive_dotState, & + constitutive_collectDotState, & + constitutive_dotTemperature implicit none @@ -258,7 +274,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) 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 + real(pReal), dimension(constitutive_maxSizeState) :: myState, & ! local copy of the state + myDotState ! local copy of dotState integer(pInt) NiterationCrystallite, & ! number of iterations in crystallite loop NiterationState ! number of iterations in state loop integer(pInt) e, & ! element index @@ -267,8 +284,11 @@ subroutine crystallite_stressAndItsTangent(updateJaco) k, & l, & myNgrains, & - mySizeState - logical onTrack, & ! flag indicating wether we are still on track + mySizeState, & + mySizeDotState + 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 @@ -313,22 +333,19 @@ subroutine crystallite_stressAndItsTangent(updateJaco) NiterationCrystallite = 0_pInt - do while (any(crystallite_subStep(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) > subStepMin)) ! cutback loop for crystallites - - NiterationCrystallite = NiterationCrystallite + 1 + do while (any(crystallite_subStep(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) > subStepMin)) ! 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 + 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)) 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_converged(g,i,e)) then if (debugger) then !$OMP CRITICAL (write2out) @@ -347,11 +364,17 @@ subroutine crystallite_stressAndItsTangent(updateJaco) crystallite_subLp0(:,:,g,i,e) = crystallite_Lp(:,:,g,i,e) ! ...plastic velocity gradient constitutive_subState0(g,i,e)%p = constitutive_state(g,i,e)%p ! ...microstructure crystallite_subTstar0_v(:,g,i,e) = crystallite_Tstar_v(:,g,i,e) ! ...2nd PK stress + else ! already at final time + !$OMP CRITICAL (distributionCrystallite) + debug_CrystalliteLoopDistribution(min(nCryst+1,NiterationCrystallite)) = & + debug_CrystalliteLoopDistribution(min(nCryst+1,NiterationCrystallite)) + 1 + !$OMPEND CRITICAL (distributionCrystallite) endif else crystallite_subStep(g,i,e) = 0.5_pReal * crystallite_subStep(g,i,e) ! cut step in half and restore... crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) ! ...temperature crystallite_Fp(:,:,g,i,e) = crystallite_subFp0(:,:,g,i,e) ! ...plastic def grad + crystallite_invFp(:,:,g,i,e) = math_inv3x3(crystallite_Fp(:,:,g,i,e)) crystallite_Lp(:,:,g,i,e) = crystallite_subLp0(:,:,g,i,e) ! ...plastic velocity grad constitutive_state(g,i,e)%p = constitutive_subState0(g,i,e)%p ! ...microstructure crystallite_Tstar_v(:,g,i,e) = crystallite_subTstar0_v(:,g,i,e) ! ...2nd PK stress @@ -367,8 +390,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) crystallite_onTrack(g,i,e) = crystallite_subStep(g,i,e) > subStepMin ! 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)) + crystallite_subStep(g,i,e) * & + (crystallite_partionedF(:,:,g,i,e) - crystallite_partionedF0(:,:,g,i,e)) crystallite_subdt(g,i,e) = crystallite_subStep(g,i,e) * crystallite_dt(g,i,e) if (debugger) then !$OMP CRITICAL (write2out) @@ -383,36 +406,54 @@ subroutine crystallite_stressAndItsTangent(updateJaco) enddo !$OMPEND PARALLEL DO + crystallite_nonfinished = ( crystallite_requested & + .and. crystallite_onTrack & + .and. .not. crystallite_converged) + ! --+>> preguess for state <<+-- ! ! incrementing by crystallite_subdt ! based on constitutive_subState0 ! results in constitutive_state + ! first loop for collection of state evolution based on old state + ! second loop for updating to new state !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed - myNgrains = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed - do g = 1,myNgrains - if ( crystallite_requested(g,i,e) & - .and. crystallite_onTrack(g,i,e) & - .and. .not. crystallite_converged(g,i,e)) then ! all undone crystallites - crystallite_converged(g,i,e) = crystallite_updateState(g,i,e) ! use former evolution rate - crystallite_converged(g,i,e) = crystallite_updateTemperature(g,i,e) ! use former evolution rate - crystallite_converged(g,i,e) = .false. ! force at least one iteration step even if state already converged - endif - enddo - enddo - enddo + 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 + constitutive_dotState(g,i,e)%p = 0.0_pReal ! zero out dotState + enddo; enddo; enddo + !$OMPEND PARALLEL DO + !$OMP PARALLEL DO + do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed + myNgrains = homogenization_Ngrains(mesh_element(3,e)) + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed + do g = 1,myNgrains + if (crystallite_nonfinished(g,i,e)) & ! all undone crystallites + call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_Fp(:,:,g,i,e), & + crystallite_invFp(:,:,g,i,e), crystallite_Temperature(g,i,e), g, i, e) + enddo; enddo; enddo + !$OMPEND PARALLEL DO + !$OMP PARALLEL DO + do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed + myNgrains = homogenization_Ngrains(mesh_element(3,e)) + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed + do g = 1,myNgrains + if (crystallite_nonfinished(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 ! --+>> state loop <<+-- NiterationState = 0_pInt - do while ( any( crystallite_requested(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) & - .and. crystallite_onTrack(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) & - .and. .not. crystallite_converged(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) ) & + do while ( any(crystallite_nonfinished(:,:,FEsolving_execELem(1):FEsolving_execElem(2))) & .and. NiterationState < nState) ! convergence loop for crystallite NiterationState = NiterationState + 1_pInt @@ -430,49 +471,69 @@ 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_requested(g,i,e) & - .and. crystallite_onTrack(g,i,e) & - .and. .not. crystallite_converged(g,i,e) ) & ! all undone crystallites + if (crystallite_nonfinished(g,i,e)) & ! all undone crystallites crystallite_onTrack(g,i,e) = crystallite_integrateStress(g,i,e) enddo enddo enddo !$OMPEND PARALLEL DO - + + crystallite_nonfinished = crystallite_nonfinished .and. crystallite_onTrack + ! --+>> state integration <<+-- ! ! incrementing by crystallite_subdt ! based on constitutive_subState0 ! results in constitutive_state + ! first loop for collection of state evolution based on old state + ! second loop for updating to new state + !$OMP PARALLEL DO + do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed + myNgrains = homogenization_Ngrains(mesh_element(3,e)) + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed + do g = 1,myNgrains + constitutive_dotState(g,i,e)%p = 0.0_pReal ! zero out dotState + enddo; enddo; enddo + !$OMPEND PARALLEL DO + !$OMP PARALLEL DO + do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed + myNgrains = homogenization_Ngrains(mesh_element(3,e)) + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed + do g = 1,myNgrains + if (crystallite_nonfinished(g,i,e)) & ! all undone crystallites + call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_Fp(:,:,g,i,e), & + crystallite_invFp(:,:,g,i,e), crystallite_Temperature(g,i,e), g, i, e) + enddo; enddo; enddo + !$OMPEND PARALLEL DO !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains debugger = (e == 1 .and. i == 1 .and. g == 1) - if ( crystallite_requested(g,i,e) & - .and. crystallite_onTrack(g,i,e) & - .and. .not. crystallite_converged(g,i,e)) then ! all undone crystallites - crystallite_converged(g,i,e) = crystallite_updateState(g,i,e) .and. & - crystallite_updateTemperature(g,i,e) + if (crystallite_nonfinished(g,i,e)) then ! all undone crystallites + crystallite_stateConverged(g,i,e) = crystallite_updateState(g,i,e) ! update state + crystallite_temperatureConverged(g,i,e) = crystallite_updateTemperature(g,i,e) ! update temperature + crystallite_converged(g,i,e) = crystallite_stateConverged(g,i,e) .and. crystallite_temperatureConverged(g,i,e) if (debugger) write (6,*) g,i,e,'converged after updState',crystallite_converged(g,i,e) if (crystallite_converged(g,i,e)) then !$OMP CRITICAL (distributionState) - debug_StateLoopDistribution(NiterationState) = debug_StateLoopDistribution(NiterationState) + 1 + debug_CrystalliteStateLoopDistribution(NiterationState) = & + debug_CrystalliteStateLoopDistribution(NiterationState) + 1 !$OMPEND CRITICAL (distributionState) - !$OMP CRITICAL (distributionCrystallite) - debug_CrystalliteLoopDistribution(NiterationCrystallite) = & - debug_CrystalliteLoopDistribution(NiterationCrystallite) + 1 - !$OMPEND CRITICAL (distributionCrystallite) endif endif enddo enddo enddo !$OMPEND PARALLEL DO - + + crystallite_nonfinished = crystallite_nonfinished .and. .not. crystallite_converged + enddo ! crystallite convergence loop + + NiterationCrystallite = NiterationCrystallite + 1 enddo ! cutback loop @@ -511,10 +572,12 @@ subroutine crystallite_stressAndItsTangent(updateJaco) ! debugger = (g == 1 .and. i == 1 .and. e == 1) 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 - myState(1:mySizeState) = constitutive_state(g,i,e)%p ! remember unperturbed, converged state... + 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) - myTemperature = crystallite_Temperature(g,i,e) myFe = crystallite_Fe(:,:,g,i,e) myLp = crystallite_Lp(:,:,g,i,e) myTstar_v = crystallite_Tstar_v(:,g,i,e) @@ -545,11 +608,14 @@ subroutine crystallite_stressAndItsTangent(updateJaco) onTrack = .true. converged = .false. NiterationState = 0_pInt - do while(.not. converged .and. onTrack .and. NiterationState < nState) ! keep cycling until done (potentially non-converged) + 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) converged = crystallite_updateState(g,i,e).AND.& ! update state - crystallite_updateTemperature(g,i,e) ! update temperature + onTrack = crystallite_integrateStress(g,i,e) ! stress of perturbed situation (overwrites _P,_Tstar_v,_Fp,_Lp,_Fe) + if (onTrack) then + 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,*) '-------------' @@ -563,16 +629,18 @@ subroutine crystallite_stressAndItsTangent(updateJaco) enddo if (converged) & ! converged state warrants stiffness update crystallite_dPdF(:,:,k,l,g,i,e) = (crystallite_P(:,:,g,i,e) - myP)/pert_Fg ! tangent dP_ij/dFg_kl - constitutive_state(g,i,e)%p = myState ! restore unperturbed, converged state... - crystallite_Temperature(g,i,e)= myTemperature ! ... temperature - crystallite_Fp(:,:,g,i,e) = myFp ! ... and kinematics - 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 + 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_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 + if (nState < NiterationState) write(6,*) 'ohh shit!! stiffenss state loop debugging exceeded',NiterationState !$OMPEND CRITICAL (out) enddo enddo @@ -633,9 +701,11 @@ endsubroutine mySize = constitutive_sizeDotState(g,i,e) ! calculate the residuum + if (any(abs(constitutive_dotState(g,i,e)%p) > 1e10_pReal)) & + write(6,*) 'dotState: crap found at',g,i,e,constitutive_dotState(g,i,e)%p 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(crystallite_Tstar_v(:,g,i,e),crystallite_Temperature(g,i,e),g,i,e) + 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 @@ -702,7 +772,8 @@ endsubroutine ! calculate the residuum call system_clock(count=tick,count_rate=tickrate,count_max=maxticks) residuum = crystallite_Temperature(g,i,e) - crystallite_subTemperature0(g,i,e) - & - crystallite_subdt(g,i,e) * constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e),crystallite_Temperature(g,i,e),g,i,e) + crystallite_subdt(g,i,e) * & + constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e),crystallite_Temperature(g,i,e),g,i,e) call system_clock(count=tock,count_rate=tickrate,count_max=maxticks) debug_cumDotTemperatureCalls = debug_cumDotTemperatureCalls + 1_pInt debug_cumDotTemperatureTicks = debug_cumDotTemperatureTicks + tock-tick @@ -844,7 +915,7 @@ endsubroutine A = math_mul33x33(transpose(invFp_current), math_mul33x33(transpose(Fg_new),math_mul33x33(Fg_new,invFp_current))) ! update microstructure - call constitutive_microstructure(crystallite_Temperature(g,i,e),g,i,e) + call constitutive_microstructure(crystallite_subTemperature0(g,i,e), crystallite_subFp0, g, i, e) ! get elasticity tensor C_66 = constitutive_homogenizedC(g,i,e) @@ -877,7 +948,7 @@ LpLoop: do ! calculate plastic velocity gradient and its tangent according to constitutive law call system_clock(count=tick,count_rate=tickrate,count_max=maxticks) - call constitutive_LpAndItsTangent(Lp_constitutive,dLp_constitutive,Tstar_v,crystallite_Temperature(g,i,e),g,i,e) + call constitutive_LpAndItsTangent(Lp_constitutive, dLp_constitutive, Tstar_v, crystallite_Temperature(g,i,e), g, i, e) call system_clock(count=tock,count_rate=tickrate,count_max=maxticks) debug_cumLpCalls = debug_cumLpCalls + 1_pInt debug_cumLpTicks = debug_cumLpTicks + tock-tick @@ -934,6 +1005,10 @@ LpLoop: do if (debugger) then !$OMP CRITICAL (write2out) write(6,*) '::: integrateStress failed on dR/dLp inversion at iteration', NiterationStress + write(6,'(a9,3(i3,x),/,9(9(e12.2,x)/))') 'dTdLp at ',g,i,e,dTdLp + write(6,'(a20,3(i3,x),/,9(9(e12.2,x)/))') 'dLp_constitutive at ',g,i,e,dLp_constitutive + write(6,'(a9,3(i3,x),/,9(9(f12.7,x)/))') 'dRdLp at ',g,i,e,dRdLp + write(6,'(a11,3(i3,x),/,3(3(f12.7,x)/))') 'Lpguess at ',g,i,e,Lpguess !$OMPEND CRITICAL (write2out) endif return @@ -979,6 +1054,7 @@ LpLoop: do crystallite_Tstar_v(:,g,i,e) = Tstar_v crystallite_Fp(:,:,g,i,e) = Fp_new crystallite_Fe(:,:,g,i,e) = Fe_new + crystallite_invFp(:,:,g,i,e) = invFp_new ! set return flag to true crystallite_integrateStress = .true. @@ -993,6 +1069,7 @@ LpLoop: do !$OMP CRITICAL (distributionStress) debug_StressLoopDistribution(NiterationStress) = debug_StressLoopDistribution(NiterationStress) + 1 + if (nStress < NiterationStress) write(6,*) 'ohh shit!! debug loop of stress exceeded',NiterationStress !$OMPEND CRITICAL (distributionStress) return diff --git a/code/debug.f90 b/code/debug.f90 index b853f9b35..652689dcf 100644 --- a/code/debug.f90 +++ b/code/debug.f90 @@ -6,10 +6,11 @@ implicit none integer(pInt), dimension(:), allocatable :: debug_StressLoopDistribution - integer(pInt), dimension(:), allocatable :: debug_StateLoopDistribution + integer(pInt), dimension(:), allocatable :: debug_CrystalliteStateLoopDistribution integer(pInt), dimension(:), allocatable :: debug_StiffnessStateLoopDistribution integer(pInt), dimension(:), allocatable :: debug_CrystalliteLoopDistribution - integer(pInt), dimension(:), allocatable :: debug_MaterialpointLoopDistribution ! added <<>> + integer(pInt), dimension(:), allocatable :: debug_MaterialpointStateLoopDistribution + integer(pInt), dimension(:), allocatable :: debug_MaterialpointLoopDistribution integer(pLongInt) :: debug_cumLpTicks = 0_pInt integer(pLongInt) :: debug_cumDotStateTicks = 0_pInt integer(pLongInt) :: debug_cumDotTemperatureTicks = 0_pInt @@ -27,18 +28,20 @@ subroutine debug_init() use numerics, only: nStress, & nState, & nCryst, & - nHomog ! added <<>> + nMPstate, & + nHomog implicit none write(6,*) write(6,*) '<<<+- debug init -+>>>' write(6,*) - allocate(debug_StressLoopDistribution(nStress)) ; debug_StressLoopDistribution = 0_pInt - allocate(debug_StateLoopDistribution(nState)) ; debug_StateLoopDistribution = 0_pInt - allocate(debug_StiffnessStateLoopDistribution(nState)) ; debug_StiffnessStateLoopDistribution = 0_pInt - allocate(debug_CrystalliteLoopDistribution(nCryst)) ; debug_CrystalliteLoopDistribution = 0_pInt - allocate(debug_MaterialpointLoopDistribution(nhomog)) ; debug_MaterialpointLoopDistribution = 0_pInt ! added <<>> + allocate(debug_StressLoopDistribution(nStress)) ; debug_StressLoopDistribution = 0_pInt + allocate(debug_CrystalliteStateLoopDistribution(nState)); debug_CrystalliteStateLoopDistribution = 0_pInt + allocate(debug_StiffnessStateLoopDistribution(nState)) ; debug_StiffnessStateLoopDistribution = 0_pInt + allocate(debug_CrystalliteLoopDistribution(nCryst+1)) ; debug_CrystalliteLoopDistribution = 0_pInt + allocate(debug_MaterialpointStateLoopDistribution(nMPstate)) ; debug_MaterialpointStateLoopDistribution = 0_pInt + allocate(debug_MaterialpointLoopDistribution(nHomog+1)) ; debug_MaterialpointLoopDistribution = 0_pInt endsubroutine !******************************************************************** @@ -49,11 +52,12 @@ subroutine debug_reset() use prec implicit none - debug_StressLoopDistribution = 0_pInt ! initialize debugging data - debug_StateLoopDistribution = 0_pInt - debug_StiffnessStateLoopDistribution = 0_pInt - debug_CrystalliteLoopDistribution = 0_pInt - debug_MaterialpointLoopDistribution = 0_pInt ! added <<>> + debug_StressLoopDistribution = 0_pInt ! initialize debugging data + debug_CrystalliteStateLoopDistribution = 0_pInt + debug_StiffnessStateLoopDistribution = 0_pInt + debug_CrystalliteLoopDistribution = 0_pInt + debug_MaterialpointStateLoopDistribution = 0_pInt + debug_MaterialpointLoopDistribution = 0_pInt debug_cumLpTicks = 0_pInt debug_cumDotStateTicks = 0_pInt debug_cumDotTemperatureTicks = 0_pInt @@ -70,9 +74,10 @@ endsubroutine use prec use numerics, only: nStress, & - nState, & - nCryst, & - nHomog ! added <<>> + nState, & + nCryst, & + nMPstate, & + nHomog implicit none integer(pInt) i,integral @@ -111,21 +116,21 @@ endsubroutine do i=1,nStress if (debug_StressLoopDistribution(i) /= 0) then integral = integral + i*debug_StressLoopDistribution(i) - write(6,'(i25,i10)') i,debug_StressLoopDistribution(i) + write(6,'(i25,x,i10)') i,debug_StressLoopDistribution(i) endif enddo - write(6,'(a15,i10,i10)') ' total',integral,sum(debug_StressLoopDistribution) + write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_StressLoopDistribution) integral = 0_pInt write(6,*) - write(6,*) 'distribution_StateLoop :' + write(6,*) 'distribution_CrystalliteStateLoop :' do i=1,nState - if (debug_StateLoopDistribution(i) /= 0) then - integral = integral + i*debug_StateLoopDistribution(i) - write(6,'(i25,i10)') i,debug_StateLoopDistribution(i) + if (debug_CrystalliteStateLoopDistribution(i) /= 0) then + integral = integral + i*debug_CrystalliteStateLoopDistribution(i) + write(6,'(i25,x,i10)') i,debug_CrystalliteStateLoopDistribution(i) endif enddo - write(6,'(a15,i10,i10)') ' total',integral,sum(debug_StateLoopDistribution) + write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_CrystalliteStateLoopDistribution) integral = 0_pInt write(6,*) @@ -133,36 +138,57 @@ endsubroutine do i=1,nState if (debug_StiffnessStateLoopDistribution(i) /= 0) then integral = integral + i*debug_StiffnessStateLoopDistribution(i) - write(6,'(i25,i10)') i,debug_StiffnessStateLoopDistribution(i) + write(6,'(i25,x,i10)') i,debug_StiffnessStateLoopDistribution(i) endif enddo - write(6,'(a15,i10,i10)') ' total',integral,sum(debug_StiffnessStateLoopDistribution) + write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_StiffnessStateLoopDistribution) integral = 0_pInt write(6,*) write(6,*) 'distribution_CrystalliteLoop :' - do i=1,nCryst + do i=1,nCryst+1 if (debug_CrystalliteLoopDistribution(i) /= 0) then integral = integral + i*debug_CrystalliteLoopDistribution(i) - write(6,'(i25,i10)') i,debug_CrystalliteLoopDistribution(i) + if (i <= nCryst) then + write(6,'(i25,x,i10)') i,debug_CrystalliteLoopDistribution(i) + else + write(6,'(i25,a1,i10)') i-1,'+',debug_CrystalliteLoopDistribution(i) + endif endif enddo - write(6,'(a15,i10,i10)') ' total',integral,sum(debug_CrystalliteLoopDistribution) + write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_CrystalliteLoopDistribution) write(6,*) !* Material point loop counter <<>> integral = 0_pInt write(6,*) - write(6,*) 'distribution_MaterialpointLoop :' - do i=1,nCryst - if (debug_MaterialpointLoopDistribution(i) /= 0) then - integral = integral + i*debug_MaterialpointLoopDistribution(i) - write(6,'(i25,i10)') i,debug_MaterialpointLoopDistribution(i) + write(6,*) 'distribution_MaterialpointStateLoop :' + do i=1,nMPstate + if (debug_MaterialpointStateLoopDistribution(i) /= 0) then + integral = integral + i*debug_MaterialpointStateLoopDistribution(i) + write(6,'(i25,x,i10)') i,debug_MaterialpointStateLoopDistribution(i) endif enddo - write(6,'(a15,i10,i10)') ' total',integral,sum(debug_MaterialpointLoopDistribution) + write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_MaterialpointStateLoopDistribution) write(6,*) + integral = 0_pInt + write(6,*) + write(6,*) 'distribution_MaterialpointLoop :' + do i=1,nHomog+1 + if (debug_MaterialpointLoopDistribution(i) /= 0) then + integral = integral + i*debug_MaterialpointLoopDistribution(i) + if (i <= nHomog) then + write(6,'(i25,x,i10)') i,debug_MaterialpointLoopDistribution(i) + else + write(6,'(i25,a1,i10)') i-1,'+',debug_MaterialpointLoopDistribution(i) + endif + endif + enddo + write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_MaterialpointLoopDistribution) + write(6,*) + + endsubroutine END MODULE debug diff --git a/code/homogenization.f90 b/code/homogenization.f90 index 42c4c951b..e9af941c8 100644 --- a/code/homogenization.f90 +++ b/code/homogenization.f90 @@ -185,9 +185,11 @@ subroutine materialpoint_stressAndItsTangent(& use prec, only: pInt, & pReal use numerics, only: subStepMin, & - nHomog + nHomog, & + nMPstate use FEsolving, only: FEsolving_execElem, & - FEsolving_execIP + FEsolving_execIP, & + terminallyIll use mesh, only: mesh_element use material, only: homogenization_Ngrains use constitutive, only: constitutive_state0, & @@ -209,14 +211,16 @@ subroutine materialpoint_stressAndItsTangent(& crystallite_partionedTstar0_v, & crystallite_dt, & crystallite_requested, & + crystallite_converged, & crystallite_stressAndItsTangent - use debug, only: debug_MaterialpointLoopdistribution ! added <<>> + use debug, only: debug_MaterialpointLoopDistribution, & + debug_MaterialpointStateLoopDistribution implicit none real(pReal), intent(in) :: dt logical, intent(in) :: updateJaco - integer(pInt) homogenization_Niteration + integer(pInt) NiterationHomog,NiterationMPstate integer(pInt) g,i,e,myNgrains ! ------ initialize to starting condition ------ @@ -255,7 +259,8 @@ subroutine materialpoint_stressAndItsTangent(& enddo !$OMP END PARALLEL DO - + NiterationHomog = 0_pInt + ! ------ cutback loop ------ do while (any(materialpoint_subStep(:,FEsolving_execELem(1):FEsolving_execElem(2)) > subStepMin)) ! cutback loop for material points @@ -285,7 +290,11 @@ subroutine materialpoint_stressAndItsTangent(& if (homogenization_sizeState(i,e) > 0_pInt) & homogenization_subState0(i,e)%p = homogenization_state(i,e)%p ! ...internal state of homog scheme materialpoint_subF0(:,:,i,e) = materialpoint_subF(:,:,i,e) ! ...def grad - + else ! already at final time + !$OMP CRITICAL (distributionHomog) + debug_MaterialpointLoopDistribution(min(nHomog+1,NiterationHomog)) = & + debug_MaterialpointLoopDistribution(min(nHomog+1,NiterationHomog)) + 1 + !$OMPEND CRITICAL (distributionHomog) endif ! materialpoint didn't converge, so we need a cutback here @@ -316,19 +325,19 @@ subroutine materialpoint_stressAndItsTangent(& !$OMP END PARALLEL DO !* Checks for cutback/substepping loops: added <<>> - write (6,'(a,/,8(L,x))') 'MP exceeds substep min',materialpoint_subStep(:,FEsolving_execELem(1):FEsolving_execElem(2)) > subStepMin - write (6,'(a,/,8(L,x))') 'MP requested',materialpoint_requested(:,FEsolving_execELem(1):FEsolving_execElem(2)) - write (6,'(a,/,8(f6.4,x))') 'MP subFrac',materialpoint_subFrac(:,FEsolving_execELem(1):FEsolving_execElem(2)) - write (6,'(a,/,8(f6.4,x))') 'MP subStep',materialpoint_subStep(:,FEsolving_execELem(1):FEsolving_execElem(2)) + ! write (6,'(a,/,8(L,x))') 'MP exceeds substep min',materialpoint_subStep(:,FEsolving_execELem(1):FEsolving_execElem(2)) > subStepMin + ! write (6,'(a,/,8(L,x))') 'MP requested',materialpoint_requested(:,FEsolving_execELem(1):FEsolving_execElem(2)) + ! write (6,'(a,/,8(f6.4,x))') 'MP subFrac',materialpoint_subFrac(:,FEsolving_execELem(1):FEsolving_execElem(2)) + ! write (6,'(a,/,8(f6.4,x))') 'MP subStep',materialpoint_subStep(:,FEsolving_execELem(1):FEsolving_execElem(2)) ! ------ convergence loop material point homogenization ------ - homogenization_Niteration = 0_pInt + NiterationMPstate = 0_pInt do while (any( materialpoint_requested(:,FEsolving_execELem(1):FEsolving_execElem(2)) & .and. .not. materialpoint_doneAndHappy(1,:,FEsolving_execELem(1):FEsolving_execElem(2)) & - ) .and. homogenization_Niteration < nHomog) ! convergence loop for materialpoint - homogenization_Niteration = homogenization_Niteration + 1 + ) .and. NiterationMPstate < nMPstate) ! convergence loop for materialpoint + NiterationMPstate = NiterationMPstate + 1 ! --+>> deformation partitioning <<+-- ! @@ -367,11 +376,15 @@ subroutine materialpoint_stressAndItsTangent(& do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed if ( materialpoint_requested(i,e) .and. & .not. materialpoint_doneAndHappy(1,i,e)) then - materialpoint_doneAndHappy(:,i,e) = homogenization_updateState(i,e) + if (.not. all(crystallite_converged(:,i,e))) then + materialpoint_doneAndHappy(:,i,e) = (/.true.,.false./) + else + materialpoint_doneAndHappy(:,i,e) = homogenization_updateState(i,e) + endif materialpoint_converged(i,e) = all(materialpoint_doneAndHappy(:,i,e)) ! converged if done and happy if (materialpoint_converged(i,e)) & ! added <<>> - debug_MaterialpointLoopdistribution(homogenization_Niteration) = & - debug_MaterialpointLoopdistribution(homogenization_Niteration) + 1 + debug_MaterialpointStateLoopdistribution(NiterationMPstate) = & + debug_MaterialpointStateLoopdistribution(NiterationMPstate) + 1 endif enddo enddo @@ -379,18 +392,26 @@ subroutine materialpoint_stressAndItsTangent(& enddo ! homogenization convergence loop + NiterationHomog = NiterationHomog +1_pInt + enddo ! cutback loop ! check for non-performer: any(.not. converged) - ! replace with elastic response ? + ! replace everybody with odd response ? + !$OMP PARALLEL DO - do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed +elementLoop: do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed - call homogenization_averageStressAndItsTangent(i,e) - call homogenization_averageTemperature(i,e) + if (materialpoint_converged(i,e)) then + call homogenization_averageStressAndItsTangent(i,e) + call homogenization_averageTemperature(i,e) + else + terminallyIll = .true. + exit elementLoop + endif enddo - enddo + enddo elementLoop !$OMP END PARALLEL DO write (6,*) 'Material Point finished' diff --git a/code/homogenization_RGC.f90 b/code/homogenization_RGC.f90 index 9f4e2f5c3..2bbd35b5d 100644 --- a/code/homogenization_RGC.f90 +++ b/code/homogenization_RGC.f90 @@ -64,8 +64,7 @@ subroutine homogenization_RGC_init(& allocate(homogenization_RGC_Ngrains(3,maxNinstance)); homogenization_RGC_Ngrains = 0_pInt allocate(homogenization_RGC_ciAlpha(3,maxNinstance)); homogenization_RGC_ciAlpha = 0.0_pReal allocate(homogenization_RGC_xiAlpha(3,maxNinstance)); homogenization_RGC_xiAlpha = 0.0_pReal - allocate(homogenization_RGC_output(maxval(homogenization_Noutput),maxNinstance)); & - homogenization_RGC_output = '' + allocate(homogenization_RGC_output(maxval(homogenization_Noutput),maxNinstance)); homogenization_RGC_output = '' rewind(file) line = '' diff --git a/code/material.config b/code/material.config index 083fcf441..393fe3741 100644 --- a/code/material.config +++ b/code/material.config @@ -6,56 +6,23 @@ type isostrain Ngrains 1 -[RGC] -type RGC -Ngrains 8 - [Taylor2] type isostrain Ngrains 2 -[Taylor4] -type isostrain -Ngrains 4 - -[Taylor10] -type isostrain -Ngrains 10 - -[Taylor100] -type isostrain -Ngrains 100 - ##################### ##################### -[TWIPSteel_Poly] -(constituent) phase 5 texture 1 fraction 1.0 +[Aluminum_Poly] +(constituent) phase 3 texture 1 fraction 1.0 -[TWIPsteel_001] -(constituent) phase 5 texture 2 fraction 1.0 +[Aluminum_001] +(constituent) phase 3 texture 2 fraction 1.0 -[TWIPsteel_011] -(constituent) phase 5 texture 3 fraction 1.0 - -[TWIPsteel_111] -(constituent) phase 5 texture 4 fraction 1.0 - -[TWIPsteel_123] -(constituent) phase 5 texture 5 fraction 1.0 - -[Alu_Polycrystal] +[Aluminum_j2] (constituent) phase 1 texture 1 fraction 1.0 -[Alu_001] -(constituent) phase 1 texture 2 fraction 1.0 - -[Alu_011] -(constituent) phase 1 texture 3 fraction 1.0 - -[Alu_111] -(constituent) phase 1 texture 4 fraction 1.0 ##################### @@ -101,11 +68,11 @@ c44 28.34e9 gdot0_slip 0.001 n_slip 20 -s0_slip 31e6 0 0 0 # per family -ssat_slip 63e6 0 0 0 # per family +tau0_slip 31e6 # per family +tausat_slip 63e6 # per family gdot0_twin 0.001 n_twin 20 -s0_twin 31e6 0 0 0 # per family +tau0_twin 31e6 # per family s_pr 0 # push-up factor for slip saturation due to twinning twin_b 0 twin_c 0 @@ -120,43 +87,27 @@ interaction_sliptwin 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 interaction_twinslip 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 interaction_twintwin 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 +[Aluminum_nonlocal] -[TWIP steel FeMnC] +constitution nonlocal +/nonlocal/ -constitution dislobased +(output) dislocation_density -(output) state_slip -(output) rateofshear_slip -(output) mfp_slip -(output) thresholdstress_slip -(output) state_twin -(output) rateofshear_twin -(output) mfp_twin -(output) thresholdstress_twin +lattice_structure fcc +Nslip 12 0 0 0 # per family -C11 175.0e9 # elastic constants in Pa -C12 115.0e9 -C44 135.0e9 -lattice_structure fcc -Nslip 12 -Ntwin 12 +c11 106.75e9 +c12 60.41e9 +c44 28.34e9 -### dislocation density-based constitutive parameters ### -burgers 2.56e-10 # Burgers vector [m] -Qedge 5.5e-19 # Activation energy for dislocation glide [J/K] (0.5*G*b^3) -grainsize 2.0e-5 # Average grain size [m] -stacksize 5.0e-8 # Twin stack mean thickness [m] -fmax 1.0 # Maximum admissible twin volume fraction -Ndot0 0.0 # Number of potential twin source per volume per time [1/m³.s] -interaction_coefficients 1.0 2.2 3.0 1.6 3.8 4.5 # Dislocation interaction coefficients -rho0 6.0e12 # Initial dislocation density [m/m³] -Cmfpslip 2.0 # Adjustable parameter controlling dislocation mean free path -Cmfptwin 1.0 # Adjustable parameter controlling twin mean free path -Cthresholdslip 0.1 # Adjustable parameter controlling threshold stress for dislocation motion -Cthresholdtwin 1.0 # Adjustable parameter controlling threshold stress for deformation twinning -Cactivolume 1.0 # Adjustable parameter controlling activation volume -Cstorage 0.02 # Adjustable parameter controlling dislocation storage -Carecovery 0.0 # Adjustable parameter controlling athermal recovery +burgers 2.56e-10 0 0 0 # Burgers vector in m +rhoEdgePos0 2.5e12 0 0 0 # Initial positive edge dislocation density in m/m**3 +rhoEdgeNeg0 2.5e12 0 0 0 # Initial negative edge dislocation density in m/m**3 +rhoScrewPos0 2.5e12 0 0 0 # Initial positive screw dislocation density in m/m**3 +rhoScrewNeg0 2.5e12 0 0 0 # Initial negative screw dislocation density in m/m**3 +v0 100 0 0 0 # initial dislocation velocity +interaction_SlipSlip 1.0 2.2 3.0 1.6 3.8 4.5 # Dislocation interaction coefficients ##################### diff --git a/code/math.f90 b/code/math.f90 index be609ff95..fb4858808 100644 --- a/code/math.f90 +++ b/code/math.f90 @@ -556,7 +556,6 @@ endsubroutine !************************************************************************** ! matrix multiplication 3x3 = 1 !************************************************************************** - PURE FUNCTION math_mul3x3(A,B) use prec, only: pReal, pInt @@ -578,7 +577,6 @@ endsubroutine !************************************************************************** ! matrix multiplication 6x6 = 1 !************************************************************************** - PURE FUNCTION math_mul6x6(A,B) use prec, only: pReal, pInt @@ -596,10 +594,10 @@ endsubroutine ENDFUNCTION + !************************************************************************** ! matrix multiplication 33x33 = 3x3 -!************************************************************************** - +!**************************************************************************´ PURE FUNCTION math_mul33x33(A,B) use prec, only: pReal, pInt @@ -635,27 +633,6 @@ endsubroutine ENDFUNCTION - -!************************************************************************** -! matrix multiplication 66x6 = 6 -!************************************************************************** - PURE FUNCTION math_mul66x6(A,B) - - use prec, only: pReal, pInt - implicit none - - integer(pInt) i - real(pReal), dimension(6,6), intent(in) :: A - real(pReal), dimension(6), intent(in) :: B - real(pReal), dimension(6) :: math_mul66x6 - - forall (i=1:6) math_mul66x6(i) = & - A(i,1)*B(1) + A(i,2)*B(2) + A(i,3)*B(3) + & - A(i,4)*B(4) + A(i,5)*B(5) + A(i,6)*B(6) - return - - ENDFUNCTION - !************************************************************************** ! matrix multiplication 99x99 = 9x9 @@ -680,6 +657,62 @@ endsubroutine ENDFUNCTION +!************************************************************************** +! matrix multiplication 33x3 = 3 +!************************************************************************** + PURE FUNCTION math_mul33x3(A,B) + + use prec, only: pReal, pInt + implicit none + + integer(pInt) i + real(pReal), dimension(3,3), intent(in) :: A + real(pReal), dimension(3), intent(in) :: B + real(pReal), dimension(3) :: math_mul33x3 + + forall (i=1:3) math_mul33x3(i) = A(i,1)*B(1) + A(i,2)*B(2) + A(i,3)*B(3) + return + + ENDFUNCTION + + +!************************************************************************** +! matrix multiplication 66x6 = 6 +!************************************************************************** + PURE FUNCTION math_mul66x6(A,B) + + use prec, only: pReal, pInt + implicit none + + integer(pInt) i + real(pReal), dimension(6,6), intent(in) :: A + real(pReal), dimension(6), intent(in) :: B + real(pReal), dimension(6) :: math_mul66x6 + + forall (i=1:6) math_mul66x6(i) = & + A(i,1)*B(1) + A(i,2)*B(2) + A(i,3)*B(3) + & + A(i,4)*B(4) + A(i,5)*B(5) + A(i,6)*B(6) + return + + ENDFUNCTION + + +!************************************************************************** +! transposition of a 3x3 matrix +!************************************************************************** +function math_transpose3x3(A) + + use prec, only: pReal,pInt + implicit none + + real(pReal),dimension(3,3),intent(in) :: A + real(pReal),dimension(3,3) :: math_transpose3x3 + + math_transpose3x3 = reshape( (/A(1,1), A(1,2), A(1,3), A(2,1), A(2,2), A(2,3), A(3,1), A(3,2), A(3,3) /),(/3,3/) ) + return + + endfunction + !************************************************************************** ! Cramer inversion of 3x3 matrix (function) @@ -1024,7 +1057,24 @@ endsubroutine ENDFUNCTION + +!******************************************************************** +! euclidic norm of a 3x1 vector +!******************************************************************** + pure function math_norm3(v3) + use prec, only: pReal,pInt + implicit none + + real(pReal), dimension(3), intent(in) :: v3 + real(pReal) math_norm3 + + math_norm3 = sqrt(v3(1)*v3(1)+v3(2)*v3(2)+v3(3)*v3(3)) + return + + endfunction + + !******************************************************************** ! convert 3x3 matrix into vector 9x1 !******************************************************************** diff --git a/code/mesh.f90 b/code/mesh.f90 index 4f168a8c0..e167b7213 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -46,14 +46,15 @@ integer(pInt), dimension(2) :: mesh_maxValStateVar = 0_pInt character(len=64), dimension(:), allocatable :: mesh_nameElemSet integer(pInt), dimension(:,:), allocatable :: mesh_mapElemSet - integer(pInt), dimension(:,:), allocatable, target :: mesh_mapFEtoCPelem,mesh_mapFEtoCPnode + integer(pInt), dimension(:,:), allocatable, target :: mesh_mapFEtoCPelem, mesh_mapFEtoCPnode integer(pInt), dimension(:,:), allocatable :: mesh_element, mesh_sharedElem integer(pInt), dimension(:,:,:,:), allocatable :: mesh_ipNeighborhood - real(pReal), dimension(:,:,:), allocatable :: mesh_subNodeCoord ! coordinates of subnodes per element - real(pReal), dimension(:,:), allocatable :: mesh_ipVolume ! volume associated with IP - real(pReal), dimension(:,:,:), allocatable :: mesh_ipArea ! area of interface to neighboring IP - real(pReal), dimension(:,:,:,:), allocatable :: mesh_ipAreaNormal ! area normal of interface to neighboring IP + real(pReal), dimension(:,:,:), allocatable :: mesh_subNodeCoord ! coordinates of subnodes per element + real(pReal), dimension(:,:), allocatable :: mesh_ipVolume ! volume associated with IP + real(pReal), dimension(:,:,:), allocatable :: mesh_ipArea, & ! area of interface to neighboring IP + mesh_ipCenterOfGravity ! center of gravity of IP + real(pReal), dimension(:,:,:,:), allocatable :: mesh_ipAreaNormal ! area normal of interface to neighboring IP real(pReal), allocatable :: mesh_node (:,:) integer(pInt), dimension(:,:,:,:), allocatable :: FE_nodesAtIP @@ -1577,7 +1578,9 @@ subroutine mesh_get_nodeElemDimensions (unit) real(pReal), dimension(Ntriangles,FE_NipFaceNodes) :: volume ! volumes of possible tetrahedra real(pReal), dimension(3) :: centerOfGravity - allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems)) ; mesh_ipVolume = 0.0_pReal + allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems)) ; mesh_ipVolume = 0.0_pReal + allocate(mesh_ipCenterOfGravity(3,mesh_maxNips,mesh_NcpElems)) ; mesh_ipCenterOfGravity = 0.0_pReal + do e = 1,mesh_NcpElems ! loop over cpElems t = mesh_element(2,e) ! get elemType do i = 1,FE_Nips(t) ! loop over IPs of elem @@ -1613,6 +1616,7 @@ subroutine mesh_get_nodeElemDimensions (unit) mesh_ipVolume(i,e) = mesh_ipVolume(i,e) + sum(volume) ! add contribution from this interface enddo mesh_ipVolume(i,e) = mesh_ipVolume(i,e) / FE_NipFaceNodes ! renormalize with interfaceNodeNum due to loop over them + mesh_ipCenterOfGravity(:,i,e) = centerOfGravity enddo enddo return @@ -1645,17 +1649,17 @@ subroutine mesh_get_nodeElemDimensions (unit) t = mesh_element(2,e) ! get elemType do i = 1,FE_Nips(t) ! loop over IPs of elem do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP - forall (n = 1:FE_NipFaceNodes) nPos(:,n) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e) + forall (n = 1:FE_NipFaceNodes) nPos(:,n) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e) forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles) ! start at each interface node and build valid triangles to cover interface normal(:,j,n) = math_vectorproduct(nPos(:,1+mod(n+j-1,FE_NipFaceNodes)) - nPos(:,n), & ! calc their normal vectors nPos(:,1+mod(n+j-0,FE_NipFaceNodes)) - nPos(:,n)) area(j,n) = dsqrt(sum(normal(:,j,n)*normal(:,j,n))) ! and area end forall - forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles, area(j,n) > 0.0_pReal) & - normal(:,j,n) = normal(:,j,n) / area(j,n) ! make unit normal - + forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles, area(j,n) > 0.0_pReal) & + normal(:,j,n) = normal(:,j,n) / area(j,n) ! make unit normal + mesh_ipArea(f,i,e) = sum(area) / (FE_NipFaceNodes*2.0_pReal) ! area of parallelograms instead of triangles - mesh_ipAreaNormal(:,f,i,e) = sum(sum(normal,3),2) / count(area > 0.0_pReal) ! average of all valid normals + mesh_ipAreaNormal(:,f,i,e) = sum(sum(normal,3),2) / count(area > 0.0_pReal) ! average of all valid normals enddo enddo enddo diff --git a/code/mpie_cpfem_marc.f90 b/code/mpie_cpfem_marc.f90 index 18982e6c2..7f85c398c 100644 --- a/code/mpie_cpfem_marc.f90 +++ b/code/mpie_cpfem_marc.f90 @@ -43,13 +43,14 @@ include "mesh.f90" ! uses prec, math, IO, FEsolving include "material.f90" ! uses prec, math, IO, mesh include "lattice.f90" ! uses prec, math, IO, material - include "constitutive_phenopowerlaw.f90" ! uses prec, math, IO, lattice, material, debug - include "constitutive_j2.f90" ! uses prec, math, IO, lattice, material, debug - include "constitutive_dislobased.f90" ! uses prec, math, IO, lattice, material, debug + include "constitutive_phenopowerlaw.f90" ! uses prec, math, IO, lattice, material, debug + include "constitutive_j2.f90" ! uses prec, math, IO, lattice, material, debug + include "constitutive_dislobased.f90" ! uses prec, math, IO, lattice, material, debug + include "constitutive_nonlocal.f90" ! uses prec, math, IO, lattice, material, debug include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug include "crystallite.f90" ! uses prec, math, IO, numerics - include "homogenization_isostrain.f90" ! uses prec, math, IO, - include "homogenization_RGC.f90" ! uses prec, math, IO, numerics, mesh: added <<>> + include "homogenization_isostrain.f90" ! uses prec, math, IO, + include "homogenization_RGC.f90" ! uses prec, math, IO, numerics, mesh: added <<>> include "homogenization.f90" ! uses prec, math, IO, numerics include "CPFEM.f90" ! uses prec, math, IO, numerics, debug, FEsolving, mesh, lattice, constitutive, crystallite @@ -166,7 +167,10 @@ subroutine hypela2(& if (inc == 0) then cycleCounter = 4 else - if (theCycle > ncycle .or. theInc /= inc) cycleCounter = 0 ! reset counter for each cutback or new inc + if (theCycle > ncycle .or. theInc /= inc) then + cycleCounter = 0 ! reset counter for each cutback or new inc + terminallyIll = .false. + endif if (theCycle /= ncycle .or. theLovl /= lovl) then cycleCounter = cycleCounter+1 ! ping pong outdatedFFN1 = .false. diff --git a/code/numerics.config b/code/numerics.config index a88bcb6a7..ccf994419 100644 --- a/code/numerics.config +++ b/code/numerics.config @@ -5,7 +5,8 @@ relevantStrain 1.0e-7 # strain increment considered signific iJacoStiffness 1 # frequency of stiffness update iJacoLpresiduum 1 # frequency of Jacobian update of residuum in Lp pert_Fg 1.0e-6 # strain perturbation for FEM Jacobi -nHomog 10 # homogenization loop limit +nHomog 20 # homogenization loop limit +nMPstate 10 # material point state loop limit nCryst 20 # crystallite loop limit (only for debugging info, real loop limit is "subStepMin") nState 10 # state loop limit nStress 40 # stress loop limit diff --git a/code/numerics.f90 b/code/numerics.f90 index 519ac7465..31e798d50 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -9,6 +9,7 @@ character(len=64), parameter :: numerics_configFile = 'numerics.config' ! name o integer(pInt) iJacoStiffness, & ! frequency of stiffness update iJacoLpresiduum, & ! frequency of Jacobian update of residuum in Lp nHomog, & ! homogenization loop limit + nMPstate, & ! materialpoint state loop limit nCryst, & ! crystallite loop limit (only for debugging info, real loop limit is "subStepMin") nState, & ! state loop limit nStress ! stress loop limit @@ -69,7 +70,8 @@ subroutine numerics_init() iJacoStiffness = 1_pInt iJacoLpresiduum = 1_pInt pert_Fg = 1.0e-6_pReal - nHomog = 10_pInt + nHomog = 20_pInt + nMPstate = 10_pInt nCryst = 20_pInt nState = 10_pInt nStress = 40_pInt @@ -111,6 +113,8 @@ subroutine numerics_init() pert_Fg = IO_floatValue(line,positions,2) case ('nhomog') nHomog = IO_intValue(line,positions,2) + case ('nmpstate') + nMPstate = IO_intValue(line,positions,2) case ('ncryst') nCryst = IO_intValue(line,positions,2) case ('nstate') @@ -160,6 +164,7 @@ subroutine numerics_init() write(6,'(a24,x,i8)') 'iJacoLpresiduum: ',iJacoLpresiduum write(6,'(a24,x,e8.1)') 'pert_Fg: ',pert_Fg write(6,'(a24,x,i8)') 'nHomog: ',nHomog + write(6,'(a24,x,i8)') 'nMPstate: ',nMPstate write(6,'(a24,x,i8)') 'nCryst: ',nCryst write(6,'(a24,x,i8)') 'nState: ',nState write(6,'(a24,x,i8)') 'nStress: ',nStress @@ -184,12 +189,13 @@ subroutine numerics_init() if (iJacoLpresiduum < 1_pInt) call IO_error(262) if (pert_Fg <= 0.0_pReal) call IO_error(263) if (nHomog < 1_pInt) call IO_error(264) + if (nMPstate < 1_pInt) call IO_error(279) !! missing in IO !! if (nCryst < 1_pInt) call IO_error(265) if (nState < 1_pInt) call IO_error(266) if (nStress < 1_pInt) call IO_error(267) if (subStepMin <= 0.0_pReal) call IO_error(268) if (rTol_crystalliteState <= 0.0_pReal) call IO_error(269) - if (rTol_crystalliteTemperature <= 0.0_pReal) call IO_error(276) + if (rTol_crystalliteTemperature <= 0.0_pReal) call IO_error(276) !! oops !! if (rTol_crystalliteStress <= 0.0_pReal) call IO_error(270) if (aTol_crystalliteStress <= 0.0_pReal) call IO_error(271) @@ -198,7 +204,7 @@ subroutine numerics_init() if (relTol_RGC <= 0.0_pReal) call IO_error(273) if (absMax_RGC <= 0.0_pReal) call IO_error(274) if (relMax_RGC <= 0.0_pReal) call IO_error(275) - if (pPert_RGC <= 0.0_pReal) call IO_error(276) + if (pPert_RGC <= 0.0_pReal) call IO_error(276) !! oops !! if (xSmoo_RGC <= 0.0_pReal) call IO_error(277) endsubroutine