diff --git a/code/CPFEM.f90 b/code/CPFEM.f90 index 1420f964c..d107bfb2a 100644 --- a/code/CPFEM.f90 +++ b/code/CPFEM.f90 @@ -286,7 +286,9 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt ! return the local stress and the jacobian from storage cauchyStress(1:ngens) = CPFEM_cs(1:ngens,IP,cp_en) jacobian(1:ngens,1:ngens) = CPFEM_dcsdE(1:ngens,1:ngens,IP,cp_en) - + ! return temperature + Temperature = materialpoint_Temperature(IP,cp_en) + return end subroutine diff --git a/code/constitutive.f90 b/code/constitutive.f90 index 47e5776c2..b95708ca9 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -29,6 +29,7 @@ CONTAINS !* - constitutive_microstructure !* - constitutive_LpAndItsTangent !* - constitutive_dotState +!* - constitutive_dotTemperature !* - constitutive_postResults !**************************************** @@ -275,6 +276,45 @@ function constitutive_dotState(Tstar_v,Temperature,ipc,ip,el) end function +function constitutive_dotTemperature(Tstar_v,Temperature,ipc,ip,el) +!********************************************************************* +!* This subroutine contains the constitutive equation for * +!* calculating the rate of change of microstructure * +!* INPUT: * +!* - Tstar_v : 2nd Piola Kirchhoff stress tensor (Mandel) * +!* - state : current microstructure * +!* - ipc : component-ID of current integration point * +!* - ip : current integration point * +!* - el : current element * +!* OUTPUT: * +!* - constitutive_dotTemperature : evolution of temperature * +!********************************************************************* + use prec, only: pReal,pInt + use material, only: phase_constitution,material_phase + use constitutive_phenomenological + use constitutive_j2 + use constitutive_dislobased + implicit none + +!* Definition of variables + integer(pInt) ipc,ip,el + real(pReal) Temperature + real(pReal), dimension(6) :: Tstar_v + real(pReal) constitutive_dotTemperature + + select case (phase_constitution(material_phase(ipc,ip,el))) + case (constitutive_phenomenological_label) + constitutive_dotTemperature = constitutive_phenomenological_dotTemperature(Tstar_v,Temperature,constitutive_state,ipc,ip,el) + case (constitutive_j2_label) + constitutive_dotTemperature = constitutive_j2_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) + + end select + return +end function + + pure function constitutive_postResults(Tstar_v,Temperature,dt,ipc,ip,el) !********************************************************************* !* return array of constitutive results * diff --git a/code/constitutive_dislobased.f90 b/code/constitutive_dislobased.f90 index 506d365c1..807bd682f 100644 --- a/code/constitutive_dislobased.f90 +++ b/code/constitutive_dislobased.f90 @@ -93,6 +93,7 @@ CONTAINS !* - constitutive_microstructure !* - constitutive_LpAndItsTangent !* - consistutive_dotState +!* - constitutive_dotTemperature !* - consistutive_postResults !**************************************** @@ -372,6 +373,7 @@ subroutine constitutive_dislobased_microstructure(Temperature,state,ipc,ip,el) real(pReal) Temperature type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state + Temperature = 298.0 matID = phase_constitutionInstance(material_phase(ipc,ip,el)) n = constitutive_dislobased_Nslip(matID) !* Quantities derived from state - slip @@ -401,10 +403,10 @@ subroutine constitutive_dislobased_microstructure(Temperature,state,ipc,ip,el) constitutive_dislobased_c3(matID)*state(ipc,ip,el)%p(4*n+i)*constitutive_dislobased_bg(matID)**2.0_pReal state(ipc,ip,el)%p(6*n+i) = & - (2.0_pReal*kB*Temperature*sqrt(state(ipc,ip,el)%p(2*n+i)))/& - (constitutive_dislobased_c1(matID)*constitutive_dislobased_c3(matID)*constitutive_dislobased_Gmod(matID)*& - state(ipc,ip,el)%p(4*n+i)*constitutive_dislobased_bg(matID)**3.0_pReal) - + (2.0_pReal*kB*Temperature)/(constitutive_dislobased_c1(matID)*constitutive_dislobased_c2(matID)*& + constitutive_dislobased_c3(matID)*constitutive_dislobased_Gmod(matID)*constitutive_dislobased_bg(matID)**3.0_pReal)*& + sqrt(state(ipc,ip,el)%p(n+i)*state(ipc,ip,el)%p(2*n+i)) + state(ipc,ip,el)%p(7*n+i) = & state(ipc,ip,el)%p(6*n+i)*constitutive_dislobased_bg(matID)*attack_frequency*state(ipc,ip,el)%p(4*n+i)*& exp(-constitutive_dislobased_Qedge(matID)/(kB*Temperature)) @@ -427,7 +429,7 @@ subroutine constitutive_dislobased_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Tempera !* - dLp_dTstar : derivative of Lp (4th-rank tensor) * !********************************************************************* use prec, only: pReal,pInt,p_vec - use math, only: math_Plain3333to99, math_mul6x6 + use math, only: math_Plain3333to99 use lattice, only: lattice_Sslip,lattice_Sslip_v use mesh, only: mesh_NcpElems,mesh_maxNips use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance @@ -446,6 +448,7 @@ subroutine constitutive_dislobased_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Tempera real(pReal), dimension(constitutive_dislobased_Nslip(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: & gdot_slip,dgdot_dtauslip,tau_slip + Temperature = 298.0 matID = phase_constitutionInstance(material_phase(ipc,ip,el)) n = constitutive_dislobased_Nslip(matID) @@ -453,10 +456,11 @@ subroutine constitutive_dislobased_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Tempera Lp = 0.0_pReal gdot_slip = 0.0_pReal do i = 1,constitutive_dislobased_Nslip(matID) - tau_slip(i) = math_mul6x6(Tstar_v,lattice_Sslip_v(:,i,constitutive_dislobased_structure(matID))) - if ((abs(tau_slip(i))-state(ipc,ip,el)%p(3*n+i))>0) & + tau_slip(i) = dot_product(Tstar_v,lattice_Sslip_v(:,i,constitutive_dislobased_structure(matID))) + + if (abs(tau_slip(i))-state(ipc,ip,el)%p(3*n+i)>0) & gdot_slip(i) = state(ipc,ip,el)%p(7*n+i)*sign(1.0_pReal,tau_slip(i))*& - sinh(((abs(tau_slip(i))-state(ipc,ip,el)%p(3*n+i))*state(ipc,ip,el)%p(5*n+i))/(kB*Temperature)) + sinh(((abs(tau_slip(i))-state(ipc,ip,el)%p(3*n+i))*state(ipc,ip,el)%p(5*n+i))/(kB*Temperature) ) Lp = Lp + gdot_slip(i)*lattice_Sslip(:,:,i,constitutive_dislobased_structure(matID)) enddo @@ -466,9 +470,11 @@ subroutine constitutive_dislobased_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Tempera dLp_dTstar = 0.0_pReal dgdot_dtauslip = 0.0_pReal do i = 1,constitutive_dislobased_Nslip(matID) + if ((abs(tau_slip(i))-state(ipc,ip,el)%p(3*n+i))>0) & dgdot_dtauslip(i) = (state(ipc,ip,el)%p(7*n+i)*state(ipc,ip,el)%p(5*n+i))/(kB*Temperature)*& - cosh(((abs(tau_slip(i))-state(ipc,ip,el)%p(3*n+i))*state(ipc,ip,el)%p(5*n+i))/(kB*Temperature)) + cosh(((abs(tau_slip(i))-state(ipc,ip,el)%p(3*n+i))*state(ipc,ip,el)%p(5*n+i))/(kB*Temperature)) + forall (k=1:3,l=1:3,m=1:3,n=1:3) & dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + & dgdot_dtauslip(i)*lattice_Sslip(k,l,i,constitutive_dislobased_structure(matID))* & @@ -506,13 +512,16 @@ function constitutive_dislobased_dotState(Tstar_v,Temperature,state,ipc,ip,el) real(pReal), dimension(constitutive_dislobased_Nslip(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: & constitutive_dislobased_dotState + Temperature = 298.0 matID = phase_constitutionInstance(material_phase(ipc,ip,el)) n = constitutive_dislobased_Nslip(matID) !* Dislocation density evolution constitutive_dislobased_dotState = 0.0_pReal do i = 1,n + tau_slip = dot_product(Tstar_v,lattice_Sslip_v(:,i,constitutive_dislobased_structure(matID))) + if (abs(tau_slip) > state(ipc,ip,el)%p(3*n+i)) then gdot_slip = state(ipc,ip,el)%p(7*n+i)*sign(1.0_pReal,tau_slip)*& sinh(((abs(tau_slip)-state(ipc,ip,el)%p(3*n+i))*state(ipc,ip,el)%p(5*n+i))/(kB*Temperature)) @@ -522,10 +531,6 @@ function constitutive_dislobased_dotState(Tstar_v,Temperature,state,ipc,ip,el) athermal_recovery = constitutive_dislobased_c7(matID)*state(ipc,ip,el)%p(i)*abs(gdot_slip) - !thermal_recovery = constitutive_dislobased_c8(matID)*abs(tau_slip)*state(ipc,ip,el)%p(i)**(2.0_pReal)*& - ! ((constitutive_dislobased_D0(matID)*constitutive_dislobased_bg(matID)**(3.0_pReal))/& - ! (kB*Temperature))*exp(-constitutive_dislobased_Qsd(matID)/(kB*Temperature)) - constitutive_dislobased_dotState(i) = locks - athermal_recovery endif enddo @@ -533,6 +538,38 @@ function constitutive_dislobased_dotState(Tstar_v,Temperature,state,ipc,ip,el) return end function + +function constitutive_dislobased_dotTemperature(Tstar_v,Temperature,state,ipc,ip,el) +!********************************************************************* +!* rate of change of microstructure * +!* INPUT: * +!* - Tstar_v : 2nd Piola Kirchhoff stress tensor (Mandel) * +!* - ipc : component-ID at current integration point * +!* - ip : current integration point * +!* - el : current element * +!* OUTPUT: * +!* - constitutive_dotTemperature : evolution of Temperature * +!********************************************************************* + use prec, only: pReal,pInt,p_vec + use lattice, only: lattice_Sslip_v + use mesh, only: mesh_NcpElems,mesh_maxNips + use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance + implicit none + +!* Definition of variables + integer(pInt) ipc,ip,el + integer(pInt) matID,i,n + real(pReal) Temperature + type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state + real(pReal), dimension(6) :: Tstar_v + real(pReal) constitutive_dislobased_dotTemperature + + constitutive_dislobased_dotTemperature = 0.0_pReal + + return +end function + + pure function constitutive_dislobased_postResults(Tstar_v,Temperature,dt,state,ipc,ip,el) !********************************************************************* !* return array of constitutive results * @@ -544,7 +581,6 @@ pure function constitutive_dislobased_postResults(Tstar_v,Temperature,dt,state,i !* - el : current element * !********************************************************************* use prec, only: pReal,pInt,p_vec - use math, only: math_mul6x6 use lattice, only: lattice_Sslip_v use mesh, only: mesh_NcpElems,mesh_maxNips use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance,phase_Noutput @@ -556,7 +592,7 @@ pure function constitutive_dislobased_postResults(Tstar_v,Temperature,dt,state,i real(pReal), dimension(6), intent(in) :: Tstar_v type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state integer(pInt) matID,o,i,c,n - real(pReal) tau_slip, active_rate + real(pReal) tau_slip real(pReal), dimension(constitutive_dislobased_sizePostResults(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: & constitutive_dislobased_postResults @@ -567,20 +603,23 @@ pure function constitutive_dislobased_postResults(Tstar_v,Temperature,dt,state,i do o = 1,phase_Noutput(material_phase(ipc,ip,el)) select case(constitutive_dislobased_output(o,matID)) + case ('dislodensity') constitutive_dislobased_postResults(c+1:c+n) = state(ipc,ip,el)%p(1:n) c = c + n + case ('rateofshear') do i = 1,n - tau_slip = math_mul6x6(Tstar_v,lattice_Sslip_v(:,i,constitutive_dislobased_structure(matID))) + tau_slip = dot_product(Tstar_v,lattice_Sslip_v(:,i,constitutive_dislobased_structure(matID))) if ((abs(tau_slip)-state(ipc,ip,el)%p(3*n+i))>0) then constitutive_dislobased_postResults(c+i) = state(ipc,ip,el)%p(7*n+i)*sign(1.0_pReal,tau_slip)*& - sinh(((abs(tau_slip)-state(ipc,ip,el)%p(3*n+i))*state(ipc,ip,el)%p(5*n+i))/(kB*Temperature)) + sinh(((abs(tau_slip)-state(ipc,ip,el)%p(3*n+i))*state(ipc,ip,el)%p(5*n+i))/(kB*298.0)) else constitutive_dislobased_postResults(c+i) = 0.0_pReal endif enddo c = c + n + end select enddo diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 681c313a7..7ae3fda49 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -25,7 +25,10 @@ real(pReal), dimension (:,:,:), allocatable :: crystallite_dt, & crystallite_subdt, & ! substepped time increment of each grain crystallite_subFrac, & ! already calculated fraction of increment crystallite_subStep, & ! size of next integration step - crystallite_Temperature ! Temp of each grain + crystallite_Temperature, & ! Temp of each grain + crystallite_Temperature0, & ! Temp of each grain at start of FE inc + crystallite_partionedTemperature0, & ! Temp of each grain at start of homog inc + crystallite_subTemperature0 ! Temp of each grain at start of crystallite inc real(pReal), dimension (:,:,:,:), allocatable :: crystallite_Tstar_v, & ! current 2nd Piola-Kirchhoff stress vector (end of converged time step) crystallite_Tstar0_v, & ! 2nd Piola-Kirchhoff stress vector at start of FE inc crystallite_partionedTstar0_v, & ! 2nd Piola-Kirchhoff stress vector at start of homog inc @@ -93,72 +96,39 @@ subroutine crystallite_init() eMax, & ! maximum number of elements myNgrains - !*** global variables ***! - ! crystallite_Fe - ! crystallite_Fp - ! crystallite_Lp - ! crystallite_F0 - ! crystallite_Fp0 - ! crystallite_Lp0 - ! crystallite_partionedF - ! crystallite_partionedF0 - ! crystallite_partionedFp0 - ! crystallite_partionedLp0 - ! crystallite_subF - ! crystallite_subF0 - ! crystallite_subFp0 - ! crystallite_subLp0 - ! crystallite_P - ! crystallite_Tstar_v - ! crystallite_Tstar0_v - ! crystallite_partionedTstar0_v - ! crystallite_subTstar0_v - ! crystallite_dPdF - ! crystallite_fallbackdPdF - ! crystallite_dt - ! crystallite_subdt - ! crystallite_subFrac - ! crystallite_subStep - ! crystallite_Temperature - ! crystallite_localConstitution - ! crystallite_requested - ! crystallite_onTrack - ! crystallite_converged - - !*** global functions or subroutines ***! - ! crystallite_stressAndItsTangent - - gMax = homogenization_maxNgrains iMax = mesh_maxNips eMax = mesh_NcpElems + allocate(crystallite_Temperature(gMax,iMax,eMax)); crystallite_Temperature = 0.0_pReal + allocate(crystallite_P(3,3,gMax,iMax,eMax)); crystallite_P = 0.0_pReal allocate(crystallite_Fe(3,3,gMax,iMax,eMax)); crystallite_Fe = 0.0_pReal allocate(crystallite_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_Temperature0(gMax,iMax,eMax)); crystallite_Temperature0 = 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_P(3,3,gMax,iMax,eMax)); crystallite_P = 0.0_pReal - allocate(crystallite_Tstar_v(6,gMax,iMax,eMax)); crystallite_Tstar_v = 0.0_pReal - allocate(crystallite_Tstar0_v(6,gMax,iMax,eMax)); crystallite_Tstar0_v = 0.0_pReal - allocate(crystallite_partionedTstar0_v(6,gMax,iMax,eMax)); crystallite_partionedTstar0_v = 0.0_pReal - allocate(crystallite_subTstar0_v(6,gMax,iMax,eMax)); crystallite_subTstar0_v = 0.0_pReal + allocate(crystallite_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_Temperature(gMax,iMax,eMax)); crystallite_Temperature = 0.0_pReal allocate(crystallite_localConstitution(gMax,iMax,eMax)); allocate(crystallite_requested(gMax,iMax,eMax)); crystallite_requested = .false. allocate(crystallite_onTrack(gMax,iMax,eMax)); crystallite_onTrack = .false. @@ -169,14 +139,16 @@ subroutine crystallite_init() myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element do g = 1,myNgrains + crystallite_partionedTemperature0(g,i,e) = crystallite_Temperature0(g,i,e) crystallite_Fp0(:,:,g,i,e) = math_EulerToR(material_EulerAngles(:,g,i,e)) ! plastic def gradient reflects init orientation crystallite_F0(:,:,g,i,e) = math_I3 crystallite_partionedFp0(:,:,g,i,e) = crystallite_Fp0(:,:,g,i,e) crystallite_partionedF0(:,:,g,i,e) = crystallite_F0(:,:,g,i,e) crystallite_partionedF(:,:,g,i,e) = crystallite_F0(:,:,g,i,e) crystallite_requested(g,i,e) = .true. - crystallite_localConstitution(g,i,e) = phase_localConstitution(material_phase(g,i,e)) - enddo + crystallite_localConstitution(g,i,e) = phase_localConstitution(material_phase(g,i,e)) + + enddo enddo enddo !$OMPEND PARALLEL DO @@ -190,17 +162,21 @@ subroutine crystallite_init() 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_Fe: ', shape(crystallite_Fe) + 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_Temperature0: ', shape(crystallite_Temperature0) + 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_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_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) @@ -215,8 +191,7 @@ subroutine crystallite_init() 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_Temperature: ', shape(crystallite_Temperature) - write(6,'(a32,x,7(i5,x))') 'crystallite_localConstitution: ', shape(crystallite_localConstitution) + 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) @@ -277,7 +252,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) !*** output variables ***! !*** local variables ***! - real(pReal), dimension(3,3) :: invFp, & ! inverse of the plastic deformation gradient + real(pReal) myTemperature ! local copy of the temperature + real(pReal), dimension(3,3) :: invFp, & ! inverse of the plastic deformation gradient Fe_guess, & ! guess for elastic deformation gradient Tstar, & ! 2nd Piola-Kirchhoff stress tensor myF, & ! local copy of the deformation gradient @@ -299,45 +275,14 @@ subroutine crystallite_stressAndItsTangent(updateJaco) logical onTrack, & ! flag indicating wether we are still on track converged ! flag indicating if iteration converged - !*** global variables ***! - ! crystallite_Fe - ! crystallite_Fp - ! crystallite_Lp - ! crystallite_partionedF - ! crystallite_partionedF0 - ! crystallite_partionedFp0 - ! crystallite_partionedLp0 - ! crystallite_subF - ! crystallite_subF0 - ! crystallite_subFp0 - ! crystallite_subLp0 - ! crystallite_P - ! crystallite_Tstar_v - ! crystallite_Tstar0_v - ! crystallite_partionedTstar0_v - ! crystallite_subTstar0_v - ! crystallite_dPdF - ! crystallite_fallbackdPdF - ! crystallite_dt - ! crystallite_subdt - ! crystallite_subFrac - ! crystallite_subStep - ! crystallite_Temperature - ! crystallite_localConstitution - ! crystallite_requested - ! crystallite_onTrack - ! crystallite_converged - - !*** global functions or subroutines ***! - ! crystallite_integrateStress - ! crystallite_updateState ! ------ initialize to starting condition ------ write (6,*) write (6,*) 'Crystallite request from Materialpoint' - write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedF0 of 1 1 1',crystallite_partionedF0(1:3,:,1,1,1) + write (6,'(a,/,(f12.7,x))') 'crystallite_partionedTemperature0 of 1 1 1',crystallite_partionedTemperature0(1,1,1) + write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedF0 of 1 1 1',crystallite_partionedF0(1:3,:,1,1,1) write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedFp0 of 1 1 1',crystallite_partionedFp0(1:3,:,1,1,1) write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedF of 1 1 1',crystallite_partionedF(1:3,:,1,1,1) write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedLp0 of 1 1 1',crystallite_partionedLp0(1:3,:,1,1,1) @@ -349,6 +294,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains if (crystallite_requested(g,i,e)) then ! initialize restoration point of ... + crystallite_subTemperature0(g,i,e) = crystallite_partionedTemperature0(g,i,e) ! ...temperature constitutive_subState0(g,i,e)%p = constitutive_partionedState0(g,i,e)%p ! ...microstructure crystallite_subFp0(:,:,g,i,e) = crystallite_partionedFp0(:,:,g,i,e) ! ...plastic def grad crystallite_subLp0(:,:,g,i,e) = crystallite_partionedLp0(:,:,g,i,e) ! ...plastic velocity grad @@ -389,7 +335,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) crystallite_subFrac(g,i,e) = crystallite_subFrac(g,i,e) + crystallite_subStep(g,i,e) crystallite_subStep(g,i,e) = min(1.0_pReal-crystallite_subFrac(g,i,e), 2.0_pReal * crystallite_subStep(g,i,e)) if (crystallite_subStep(g,i,e) > subStepMin) then - crystallite_subF0(:,:,g,i,e) = crystallite_subF(:,:,g,i,e) ! wind forward... + crystallite_subTemperature0(g,i,e) = crystallite_Temperature(g,i,e) ! wind forward... + crystallite_subF0(:,:,g,i,e) = crystallite_subF(:,:,g,i,e) ! ...def grad crystallite_subFp0(:,:,g,i,e) = crystallite_Fp(:,:,g,i,e) ! ...plastic def grad crystallite_subLp0(:,:,g,i,e) = crystallite_Lp(:,:,g,i,e) ! ...plastic velocity gradient constitutive_subState0(g,i,e)%p = constitutive_state(g,i,e)%p ! ...microstructure @@ -405,7 +352,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) endif else crystallite_subStep(g,i,e) = 0.5_pReal * crystallite_subStep(g,i,e) ! cut step in half and restore... - crystallite_Fp(:,:,g,i,e) = crystallite_subFp0(:,:,g,i,e) ! ...plastic def grad + 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_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 @@ -453,7 +401,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) 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) + 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 @@ -507,7 +456,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) 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) + crystallite_converged(g,i,e) = crystallite_updateState(g,i,e).AND.& + crystallite_updateTemperature(g,i,e) if (crystallite_converged(g,i,e)) then !$OMP CRITICAL (distributionState) debug_StateLoopDistribution(NiterationState) = debug_StateLoopDistribution(NiterationState) + 1 @@ -567,7 +517,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) 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... myF = crystallite_subF(:,:,g,i,e) ! ... and kinematics - myFp = crystallite_Fp(:,:,g,i,e) + 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) @@ -598,7 +549,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) NiterationState = NiterationState + 1_pInt if (debugger) write (6,'(a4,x,i6)') 'loop',NiterationState 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) ! update state + if (onTrack) converged = crystallite_updateState(g,i,e).AND.& ! update state + crystallite_updateTemperature(g,i,e) ! update temperature if (debugger) then write (6,*) '-------------' write (6,'(l,x,l)') onTrack,converged @@ -611,7 +563,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) 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_Fp(:,:,g,i,e) = myFp ! ... and kinematics + 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 @@ -640,82 +593,147 @@ endsubroutine -!******************************************************************** -! update the internal state of the constitutive law -! and tell whether state has converged -!******************************************************************** - function crystallite_updateState(& - g,& ! grain number - i,& ! integration point number - e & ! element number - ) - - !*** variables and functions from other modules ***! - use prec, only: pReal, & - pInt, & - pLongInt - use numerics, only: rTol_crystalliteState - use constitutive, only: constitutive_dotState, & - constitutive_sizeDotState, & - constitutive_subState0, & - constitutive_state - use debug, only: debugger, & - debug_cumDotStateCalls, & - debug_cumDotStateTicks - - !*** input variables ***! - integer(pInt), intent(in):: e, & ! element index - i, & ! integration point index - g ! grain index - - !*** output variables ***! - logical crystallite_updateState ! flag indicating if integration suceeded - - !*** local variables ***! - real(pReal), dimension(constitutive_sizeDotState(g,i,e)) :: residuum ! residuum from evolution of microstructure - integer(pInt) mySize - integer(pLongInt) tick, & - tock, & - tickrate, & - maxticks - - !*** global variables ***! - ! crystallite_Tstar_v - ! crystallite_subdt - ! crystallite_Temperature - - mySize = constitutive_sizeDotState(g,i,e) - - ! calculate the residuum - call system_clock(count=tick,count_rate=tickrate,count_max=maxticks) - residuum = constitutive_state(g,i,e)%p(1:mySize) - constitutive_subState0(g,i,e)%p(1:mySize) - & - crystallite_subdt(g,i,e) * constitutive_dotState(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_cumDotStateCalls = debug_cumDotStateCalls + 1_pInt - debug_cumDotStateTicks = debug_cumDotStateTicks + tock-tick - if (tock < tick) debug_cumDotStateTicks = debug_cumDotStateTicks + maxticks - - ! if NaN occured then return without changing the state +!******************************************************************** +! update the internal state of the constitutive law +! and tell whether state has converged +!******************************************************************** + function crystallite_updateState(& + g,& ! grain number + i,& ! integration point number + e & ! element number + ) + + !*** variables and functions from other modules ***! + use prec, only: pReal, & + pInt, & + pLongInt + use numerics, only: rTol_crystalliteState + use constitutive, only: constitutive_dotState, & + constitutive_sizeDotState, & + constitutive_subState0, & + constitutive_state + use debug, only: debugger, & + debug_cumDotStateCalls, & + debug_cumDotStateTicks + + !*** input variables ***! + integer(pInt), intent(in):: e, & ! element index + i, & ! integration point index + g ! grain index + + !*** output variables ***! + logical crystallite_updateState ! flag indicating if integration suceeded + + !*** local variables ***! + real(pReal), dimension(constitutive_sizeDotState(g,i,e)) :: residuum ! residuum from evolution of microstructure + integer(pInt) mySize + integer(pLongInt) tick, & + tock, & + tickrate, & + maxticks + + !*** global variables ***! + ! crystallite_Tstar_v + ! crystallite_subdt + ! crystallite_Temperature + + mySize = constitutive_sizeDotState(g,i,e) + + ! calculate the residuum + call system_clock(count=tick,count_rate=tickrate,count_max=maxticks) + residuum = constitutive_state(g,i,e)%p(1:mySize) - constitutive_subState0(g,i,e)%p(1:mySize) - & + crystallite_subdt(g,i,e) * constitutive_dotState(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_cumDotStateCalls = debug_cumDotStateCalls + 1_pInt + debug_cumDotStateTicks = debug_cumDotStateTicks + tock-tick + if (tock < tick) debug_cumDotStateTicks = debug_cumDotStateTicks + maxticks + + ! if NaN occured then return without changing the state if (any(residuum/=residuum)) then - crystallite_updateState = .false. ! indicate state update failed - if (debugger) write(6,*) '::: updateState encountered NaN' - return - endif - - ! update the microstructure - constitutive_state(g,i,e)%p(1:mySize) = constitutive_state(g,i,e)%p(1:mySize) - residuum - - ! setting flag to true if state is below relative Tolerance, otherwise set it to false - crystallite_updateState = maxval(abs(residuum/constitutive_state(g,i,e)%p(1:mySize)), & - constitutive_state(g,i,e)%p(1:mySize) /= 0.0_pReal) < rTol_crystalliteState - - if (debugger) write(6,'(a,/,f12.4)') 'updated state: ', constitutive_state(g,i,e)%p(1) - - return - - endfunction - - + crystallite_updateState = .false. ! indicate state update failed + if (debugger) write(6,*) '::: updateState encountered NaN' + return + endif + + ! update the microstructure + constitutive_state(g,i,e)%p(1:mySize) = constitutive_state(g,i,e)%p(1:mySize) - residuum + + ! setting flag to true if state is below relative Tolerance, otherwise set it to false + crystallite_updateState = maxval(abs(residuum/constitutive_state(g,i,e)%p(1:mySize)), & + constitutive_state(g,i,e)%p(1:mySize) /= 0.0_pReal) < rTol_crystalliteState + + if (debugger) write(6,'(a,/,f12.4)') 'updated state: ', constitutive_state(g,i,e)%p(1) + + return + + endfunction + + +!******************************************************************** +! update the temperature of the grain +! and tell whether it has converged +!******************************************************************** + function crystallite_updateTemperature(& + g,& ! grain number + i,& ! integration point number + e & ! element number + ) + + !*** variables and functions from other modules ***! + use prec, only: pReal, & + pInt, & + pLongInt + use numerics, only: rTol_crystalliteTemperature + use constitutive, only: constitutive_dotTemperature + use debug, only: debugger, & + debug_cumDotTemperatureCalls, & + debug_cumDotTemperatureTicks + + !*** input variables ***! + integer(pInt), intent(in):: e, & ! element index + i, & ! integration point index + g ! grain index + + !*** output variables ***! + logical crystallite_updateTemperature ! flag indicating if integration suceeded + + !*** local variables ***! + real(pReal) residuum ! residuum from evolution of temperature + integer(pLongInt) tick, & + tock, & + tickrate, & + maxticks + + ! calculate the residuum + call system_clock(count=tick,count_rate=tickrate,count_max=maxticks) + residuum = crystallite_Temperature(g,i,e) - crystallite_subTemperature0(g,i,e) - & + crystallite_subdt(g,i,e) * constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e),crystallite_Temperature(g,i,e),g,i,e) + call system_clock(count=tock,count_rate=tickrate,count_max=maxticks) + debug_cumDotTemperatureCalls = debug_cumDotTemperatureCalls + 1_pInt + debug_cumDotTemperatureTicks = debug_cumDotTemperatureTicks + tock-tick + if (tock < tick) debug_cumDotTemperatureTicks = debug_cumDotTemperatureTicks + maxticks + + ! if NaN occured then return without changing the state + if (residuum/=residuum) then + crystallite_updateTemperature = .false. ! indicate update failed + if (debugger) write(6,*) '::: updateTemperature encountered NaN' + return + endif + + ! update the microstructure + crystallite_Temperature(g,i,e) = crystallite_Temperature(g,i,e) - residuum + + ! setting flag to true if state is below relative Tolerance, otherwise set it to false + crystallite_updateTemperature = maxval(abs(residuum/crystallite_Temperature(g,i,e)), & + crystallite_Temperature(g,i,e) /= 0.0_pReal) < rTol_crystalliteTemperature + + if (debugger) write(6,'(a,/,f12.4)') 'updated temperature: ', crystallite_Temperature(g,i,e) + + return + + endfunction + + !*********************************************************************** !*** calculation of stress (P) with time integration *** diff --git a/code/debug.f90 b/code/debug.f90 index c10b3de61..95e71e507 100644 --- a/code/debug.f90 +++ b/code/debug.f90 @@ -11,8 +11,10 @@ integer(pInt), dimension(:), allocatable :: debug_CrystalliteLoopDistribution integer(pLongInt) :: debug_cumLpTicks = 0_pInt integer(pLongInt) :: debug_cumDotStateTicks = 0_pInt + integer(pLongInt) :: debug_cumDotTemperatureTicks = 0_pInt integer(pInt) :: debug_cumLpCalls = 0_pInt integer(pInt) :: debug_cumDotStateCalls = 0_pInt + integer(pInt) :: debug_cumDotTemperatureCalls = 0_pInt logical :: debugger = .false. logical :: distribution_init = .false. @@ -50,8 +52,10 @@ subroutine debug_reset() debug_CrystalliteLoopDistribution = 0_pInt debug_cumLpTicks = 0_pInt debug_cumDotStateTicks = 0_pInt + debug_cumDotTemperatureTicks = 0_pInt debug_cumLpCalls = 0_pInt debug_cumDotStateCalls = 0_pInt + debug_cumDotTemperatureCalls = 0_pInt endsubroutine @@ -87,6 +91,14 @@ endsubroutine dble(debug_cumDotStateTicks)/tickrate/1.0e-6_pReal/debug_cumDotStateCalls write(6,'(a33,x,i12)') 'total CPU ticks :',debug_cumDotStateTicks endif + write(6,*) + write(6,'(a33,x,i12)') 'total calls to dotTemperature :',debug_cumDotTemperatureCalls + if (debug_cumdotTemperatureCalls > 0_pInt) then + call system_clock(count_rate=tickrate) + write(6,'(a33,x,f12.6)') 'avg CPU time/microsecs per call :',& + dble(debug_cumDotTemperatureTicks)/tickrate/1.0e-6_pReal/debug_cumDotTemperatureCalls + write(6,'(a33,x,i12)') 'total CPU ticks :',debug_cumDotTemperatureTicks + endif integral = 0_pInt write(6,*) diff --git a/code/homogenization.f90 b/code/homogenization.f90 index aaf3301fc..d706ae69b 100644 --- a/code/homogenization.f90 +++ b/code/homogenization.f90 @@ -180,13 +180,16 @@ subroutine materialpoint_stressAndItsTangent(& use constitutive, only: constitutive_state0, & constitutive_partionedState0, & constitutive_state - use crystallite, only: crystallite_F0, & + use crystallite, only: crystallite_Temperature0, & + crystallite_Temperature, & + crystallite_F0, & crystallite_Fp0, & crystallite_Fp, & crystallite_Lp0, & crystallite_Lp, & crystallite_Tstar0_v, & - crystallite_Tstar_v, & + crystallite_Tstar_v, & + crystallite_partionedTemperature0, & crystallite_partionedF0, & crystallite_partionedF, & crystallite_partionedFp0, & @@ -218,11 +221,12 @@ subroutine materialpoint_stressAndItsTangent(& do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed ! initialize restoration points of grain... - forall (g = 1:myNgrains) constitutive_partionedState0(g,i,e)%p = constitutive_state0(g,i,e)%p ! ...microstructures - crystallite_partionedFp0(:,:,1:myNgrains,i,e) = crystallite_Fp0(:,:,1:myNgrains,i,e) ! ...plastic def grads - crystallite_partionedLp0(:,:,1:myNgrains,i,e) = crystallite_Lp0(:,:,1:myNgrains,i,e) ! ...plastic velocity grads - crystallite_partionedF0(:,:,1:myNgrains,i,e) = crystallite_F0(:,:,1:myNgrains,i,e) ! ...def grads - crystallite_partionedTstar0_v(:,1:myNgrains,i,e)= crystallite_Tstar0_v(:,1:myNgrains,i,e) ! ...2nd PK stress + forall (g = 1:myNgrains) constitutive_partionedState0(g,i,e)%p = constitutive_state0(g,i,e)%p ! ...microstructures + crystallite_partionedTemperature0(1:myNgrains,i,e) = crystallite_Temperature0(1:myNgrains,i,e)! ...temperatures + crystallite_partionedFp0(:,:,1:myNgrains,i,e) = crystallite_Fp0(:,:,1:myNgrains,i,e) ! ...plastic def grads + crystallite_partionedLp0(:,:,1:myNgrains,i,e) = crystallite_Lp0(:,:,1:myNgrains,i,e) ! ...plastic velocity grads + crystallite_partionedF0(:,:,1:myNgrains,i,e) = crystallite_F0(:,:,1:myNgrains,i,e) ! ...def grads + crystallite_partionedTstar0_v(:,1:myNgrains,i,e)= crystallite_Tstar0_v(:,1:myNgrains,i,e) ! ...2nd PK stress ! initialize restoration points of ... if (homogenization_sizeState(i,e) > 0_pInt) & @@ -258,7 +262,8 @@ subroutine materialpoint_stressAndItsTangent(& if (materialpoint_subStep(i,e) > subStepMin) then ! wind forward grain starting point of... - crystallite_partionedF0(:,:,1:myNgrains,i,e) = crystallite_partionedF(:,:,1:myNgrains,i,e) ! ...def grads + crystallite_partionedTemperature0(1:myNgrains,i,e) = crystallite_Temperature(1:myNgrains,i,e) ! ...temperatures + crystallite_partionedF0(:,:,1:myNgrains,i,e) = crystallite_partionedF(:,:,1:myNgrains,i,e) ! ...def grads crystallite_partionedFp0(:,:,1:myNgrains,i,e) = crystallite_Fp(:,:,1:myNgrains,i,e) ! ...plastic def grads crystallite_partionedLp0(:,:,1:myNgrains,i,e) = crystallite_Lp(:,:,1:myNgrains,i,e) ! ...plastic velocity grads crystallite_partionedTstar0_v(:,1:myNgrains,i,e) = crystallite_Tstar_v(:,1:myNgrains,i,e) ! ...2nd PK stress @@ -275,7 +280,8 @@ subroutine materialpoint_stressAndItsTangent(& materialpoint_subStep(i,e) = 0.5_pReal * materialpoint_subStep(i,e) ! restore... - crystallite_Fp(:,:,1:myNgrains,i,e) = crystallite_partionedFp0(:,:,1:myNgrains,i,e) ! ...plastic def grads + crystallite_Temperature(1:myNgrains,i,e) = crystallite_partionedTemperature0(1:myNgrains,i,e) ! ...temperatures + crystallite_Fp(:,:,1:myNgrains,i,e) = crystallite_partionedFp0(:,:,1:myNgrains,i,e) ! ...plastic def grads crystallite_Lp(:,:,1:myNgrains,i,e) = crystallite_partionedLp0(:,:,1:myNgrains,i,e) ! ...plastic velocity grads crystallite_Tstar_v(:,1:myNgrains,i,e) = crystallite_partionedTstar0_v(:,1:myNgrains,i,e) ! ...2nd PK stress forall (g = 1:myNgrains) constitutive_state(g,i,e)%p = constitutive_partionedState0(g,i,e)%p ! ...microstructures @@ -341,7 +347,7 @@ 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) + materialpoint_doneAndHappy(:,i,e) = homogenization_updateState(i,e) materialpoint_converged(i,e) = all(materialpoint_doneAndHappy(:,i,e)) ! converged if done and happy endif enddo @@ -358,7 +364,8 @@ subroutine materialpoint_stressAndItsTangent(& !$OMP PARALLEL DO 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_averageStressAndItsTangent(i,e) + call homogenization_averageTemperature(i,e) enddo enddo !$OMP END PARALLEL DO @@ -476,38 +483,65 @@ function homogenization_updateState(& endfunction -!******************************************************************** -! derive average stress and stiffness from constituent quantities -!******************************************************************** -subroutine homogenization_averageStressAndItsTangent(& - ip, & ! integration point - el & ! element - ) - use prec, only: pReal,pInt - use mesh, only: mesh_element - use material, only: homogenization_type, homogenization_maxNgrains - use crystallite, only: crystallite_P,crystallite_dPdF - - use homogenization_isostrain - implicit none - - integer(pInt), intent(in) :: ip,el - - select case(homogenization_type(mesh_element(3,el))) - case (homogenization_isostrain_label) - call homogenization_isostrain_averageStressAndItsTangent( materialpoint_P(:,:,ip,el), & - materialpoint_dPdF(:,:,:,:,ip,el),& - crystallite_P(:,:,:,ip,el), & - crystallite_dPdF(:,:,:,:,:,ip,el), & - ip, & - el) - end select - - return - -endsubroutine - - +!******************************************************************** +! derive average stress and stiffness from constituent quantities +!******************************************************************** +subroutine homogenization_averageStressAndItsTangent(& + ip, & ! integration point + el & ! element + ) + use prec, only: pReal,pInt + use mesh, only: mesh_element + use material, only: homogenization_type, homogenization_maxNgrains + use crystallite, only: crystallite_P,crystallite_dPdF + + use homogenization_isostrain + implicit none + + integer(pInt), intent(in) :: ip,el + + select case(homogenization_type(mesh_element(3,el))) + case (homogenization_isostrain_label) + call homogenization_isostrain_averageStressAndItsTangent( materialpoint_P(:,:,ip,el), & + materialpoint_dPdF(:,:,:,:,ip,el),& + crystallite_P(:,:,:,ip,el), & + crystallite_dPdF(:,:,:,:,:,ip,el), & + ip, & + el) + end select + + return + +endsubroutine + + +!******************************************************************** +! derive average stress and stiffness from constituent quantities +!******************************************************************** +subroutine homogenization_averageTemperature(& + ip, & ! integration point + el & ! element + ) + use prec, only: pReal,pInt + use mesh, only: mesh_element + use material, only: homogenization_type, homogenization_maxNgrains + use crystallite, only: crystallite_Temperature + + use homogenization_isostrain + implicit none + + integer(pInt), intent(in) :: ip,el + + select case(homogenization_type(mesh_element(3,el))) + case (homogenization_isostrain_label) + materialpoint_Temperature(ip,el) = homogenization_isostrain_averageTemperature(crystallite_Temperature(:,ip,el), ip, el) + end select + + return + +endsubroutine + + !******************************************************************** ! return array of homogenization results for post file inclusion ! call only, if homogenization_sizePostResults(ip,el) > 0 !! diff --git a/code/homogenization_isostrain.f90 b/code/homogenization_isostrain.f90 index fdae1ad97..64eea767d 100644 --- a/code/homogenization_isostrain.f90 +++ b/code/homogenization_isostrain.f90 @@ -232,7 +232,36 @@ subroutine homogenization_isostrain_averageStressAndItsTangent(& return -endsubroutine +endsubroutine + + +!******************************************************************** +! derive average stress and stiffness from constituent quantities +!******************************************************************** +function homogenization_isostrain_averageTemperature(& + Temperature, & ! temperature + ip, & ! my integration point + el & ! my element + ) + + use prec, only: pReal,pInt,p_vec + use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips + use material, only: homogenization_maxNgrains, homogenization_Ngrains + implicit none + +!* Definition of variables + real(pReal), dimension (homogenization_maxNgrains), intent(in) :: Temperature + integer(pInt), intent(in) :: ip,el + real(pReal) homogenization_isostrain_averageTemperature + integer(pInt) homID, i, Ngrains + +! homID = homogenization_typeInstance(mesh_element(3,el)) + Ngrains = homogenization_Ngrains(mesh_element(3,el)) + homogenization_isostrain_averageTemperature = sum(Temperature)/Ngrains + + return + +endfunction !********************************************************************