diff --git a/code/CPFEM.f90 b/code/CPFEM.f90 index d107bfb2a..9240e02fd 100644 --- a/code/CPFEM.f90 +++ b/code/CPFEM.f90 @@ -121,8 +121,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt integer(pInt), intent(in) :: element, & ! FE element number IP, & ! FE integration point number ngens ! size of stress strain law - real(pReal), intent(in) :: Temperature, & ! temperature - dt ! time increment + real(pReal), intent(inout) :: Temperature ! temperature + real(pReal), intent(in) :: dt ! time increment real(pReal), dimension (3,3), intent(in) :: ffn, & ! deformation gradient for t=t0 ffn1 ! deformation gradient for t=t1 integer(pInt), intent(in) :: mode ! computation mode 1: regular computation with aged results @@ -170,8 +170,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt call lattice_init() call material_init() call constitutive_init() - call crystallite_init() - call homogenization_init() + call crystallite_init(Temperature) ! (have to) use temperature of first IP for whole model + call homogenization_init(Temperature) call CPFEM_init() CPFEM_init_done = .true. endif @@ -262,6 +262,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt ! --+>> COLLECTION OF FEM DATA AND RETURN OF ODD STRESS AND JACOBIAN <<+-- case (3) + if (IP==1.AND.cp_en==1) write(6,*) 'Temp from CPFEM', Temperature materialpoint_Temperature(IP,cp_en) = Temperature materialpoint_F0(:,:,IP,cp_en) = ffn materialpoint_F(:,:,IP,cp_en) = ffn1 @@ -287,7 +288,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt 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) + if (theInc > 0_pInt) Temperature = materialpoint_Temperature(IP,cp_en) ! homogenized result except for potentially non-isothermal starting condition. return diff --git a/code/constitutive_dislobased.f90 b/code/constitutive_dislobased.f90 index 807bd682f..7d24b6017 100644 --- a/code/constitutive_dislobased.f90 +++ b/code/constitutive_dislobased.f90 @@ -373,7 +373,6 @@ 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 @@ -448,7 +447,6 @@ 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) @@ -512,7 +510,6 @@ 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) @@ -605,7 +602,7 @@ pure function constitutive_dislobased_postResults(Tstar_v,Temperature,dt,state,i select case(constitutive_dislobased_output(o,matID)) case ('dislodensity') - constitutive_dislobased_postResults(c+1:c+n) = state(ipc,ip,el)%p(1:n) + constitutive_dislobased_postResults(c+1:c+n) = state(ipc,ip,el)%p(6*n+1:7*n) c = c + n case ('rateofshear') @@ -613,7 +610,7 @@ pure function constitutive_dislobased_postResults(Tstar_v,Temperature,dt,state,i 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*298.0)) + sinh(((abs(tau_slip)-state(ipc,ip,el)%p(3*n+i))*state(ipc,ip,el)%p(5*n+i))/(kB*Temperature)) else constitutive_dislobased_postResults(c+i) = 0.0_pReal endif diff --git a/code/constitutive_j2.f90 b/code/constitutive_j2.f90 index b745bb431..06550c7e2 100644 --- a/code/constitutive_j2.f90 +++ b/code/constitutive_j2.f90 @@ -414,6 +414,36 @@ pure function constitutive_j2_dotState(Tstar_v, Temperature, state, g, ip, el) return +endfunction + + +!**************************************************************** +!* calculates the rate of change of microstructure * +!**************************************************************** +pure function constitutive_j2_dotTemperature(Tstar_v, Temperature, state, g, ip, el) + + !*** variables and functions from other modules ***! + 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 + + !*** input variables ***! + real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola Kirchhoff stress tensor in Mandel notation + real(pReal), intent(in) :: Temperature + integer(pInt), intent(in):: g, & ! grain number + ip, & ! integration point number + el ! element number + type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state ! state of the current microstructure + + !*** output variables ***! + real(pReal) constitutive_j2_dotTemperature ! evolution of state variable + + ! calculate dotTemperature + constitutive_j2_dotTemperature = 0.0_pReal + + return endfunction diff --git a/code/constitutive_phenomenological.f90 b/code/constitutive_phenomenological.f90 index 28c67fc92..7fe2f6581 100644 --- a/code/constitutive_phenomenological.f90 +++ b/code/constitutive_phenomenological.f90 @@ -354,19 +354,23 @@ subroutine constitutive_phenomenological_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,T !* Calculation of Lp Lp = 0.0_pReal do i = 1,constitutive_phenomenological_Nslip(matID) - tau_slip(i) = dot_product(Tstar_v,lattice_Sslip_v(:,i,constitutive_phenomenological_structure(matID))) + tau_slip(i) = dot_product(Tstar_v,lattice_Sslip_v(:,i,constitutive_phenomenological_structure(matID))) + gdot_slip(i) = constitutive_phenomenological_gdot0_slip(matID)*(abs(tau_slip(i))/state(ipc,ip,el)%p(i))**& constitutive_phenomenological_n_slip(matID)*sign(1.0_pReal,tau_slip(i)) + Lp = Lp + gdot_slip(i)*lattice_Sslip(:,:,i,constitutive_phenomenological_structure(matID)) enddo !* Calculation of the tangent of Lp dLp_dTstar3333 = 0.0_pReal dLp_dTstar = 0.0_pReal - do i = 1,constitutive_phenomenological_Nslip(matID) + do i = 1,constitutive_phenomenological_Nslip(matID) + dgdot_dtauslip(i) = constitutive_phenomenological_gdot0_slip(matID)*(abs(tau_slip(i))/state(ipc,ip,el)%p(i))**& (constitutive_phenomenological_n_slip(matID)-1.0_pReal)*& - constitutive_phenomenological_n_slip(matID)/state(ipc,ip,el)%p(i) + constitutive_phenomenological_n_slip(matID)/state(ipc,ip,el)%p(i) + 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_phenomenological_structure(matID))* & @@ -408,10 +412,13 @@ function constitutive_phenomenological_dotState(Tstar_v,Temperature,state,ipc,ip n = constitutive_phenomenological_Nslip(matID) !* Self-Hardening of each system - do i = 1,n - tau_slip = dot_product(Tstar_v,lattice_Sslip_v(:,i,constitutive_phenomenological_structure(matID))) + do i = 1,n + + tau_slip = dot_product(Tstar_v,lattice_Sslip_v(:,i,constitutive_phenomenological_structure(matID))) + gdot_slip = constitutive_phenomenological_gdot0_slip(matID)*(abs(tau_slip)/state(ipc,ip,el)%p(i))**& - constitutive_phenomenological_n_slip(matID)*sign(1.0_pReal,tau_slip) + constitutive_phenomenological_n_slip(matID)*sign(1.0_pReal,tau_slip) + self_hardening(i) = constitutive_phenomenological_h0(matID)*(1.0_pReal-state(ipc,ip,el)%p(i)/& constitutive_phenomenological_s_sat(matID))**constitutive_phenomenological_w0(matID)*abs(gdot_slip) enddo @@ -423,6 +430,37 @@ function constitutive_phenomenological_dotState(Tstar_v,Temperature,state,ipc,ip end function + +function constitutive_phenomenological_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_phenomenological_dotTemperature + + constitutive_phenomenological_dotTemperature = 0.0_pReal + + return +end function + pure function constitutive_phenomenological_postResults(Tstar_v,Temperature,dt,state,ipc,ip,el) !********************************************************************* @@ -456,10 +494,12 @@ pure function constitutive_phenomenological_postResults(Tstar_v,Temperature,dt,s constitutive_phenomenological_postResults = 0.0_pReal do o = 1,phase_Noutput(material_phase(ipc,ip,el)) - select case(constitutive_phenomenological_output(o,matID)) + select case(constitutive_phenomenological_output(o,matID)) + case ('slipresistance') constitutive_phenomenological_postResults(c+1:c+n) = state(ipc,ip,el)%p(1:n) - c = c + n + c = c + n + case ('rateofshear') do i = 1,n tau_slip = dot_product(Tstar_v,lattice_Sslip_v(:,i,constitutive_phenomenological_structure(matID))) @@ -467,7 +507,8 @@ pure function constitutive_phenomenological_postResults(Tstar_v,Temperature,dt,s (abs(tau_slip)/state(ipc,ip,el)%p(i))**& constitutive_phenomenological_n_slip(matID) enddo - c = c + n + c = c + n + end select enddo diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 7ae3fda49..57474e901 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -19,42 +19,41 @@ implicit none ! **************************************************************** ! *** General variables for the crystallite calculation *** ! **************************************************************** -integer(pInt), parameter :: crystallite_Nresults = 5_pInt ! phaseID, volume, Euler angles +integer(pInt), parameter :: crystallite_Nresults = 5_pInt ! phaseID, volume, Euler angles -real(pReal), dimension (:,:,:), allocatable :: crystallite_dt, & ! requested time increment of each grain - 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_Temperature0, & ! Temp of each grain at start of FE inc +real(pReal), dimension (:,:,:), allocatable :: crystallite_dt, & ! requested time increment of each grain + 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_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 - 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_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 - crystallite_F0, & ! def grad at start of FE inc - crystallite_partionedF, & ! def grad to be reached at end of homog inc - crystallite_partionedF0, & ! def grad at start of homog inc - crystallite_subF, & ! def grad to be reached at end of crystallite inc - crystallite_subF0, & ! def grad at start of crystallite inc - crystallite_Lp, & ! current plastic velocitiy grad (end of converged time step) - crystallite_Lp0, & ! plastic velocitiy grad at start of FE inc - crystallite_partionedLp0,& ! plastic velocity grad at start of homog inc - crystallite_subLp0,& ! plastic velocity grad at start of crystallite inc - crystallite_P ! 1st Piola-Kirchhoff stress per grain -real(pReal), dimension (:,:,:,:,:,:,:), allocatable :: crystallite_dPdF, & ! individual dPdF per grain - crystallite_fallbackdPdF ! dPdF fallback for non-converged grains (elastic prediction) + 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 + 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_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 + crystallite_F0, & ! def grad at start of FE inc + crystallite_partionedF, & ! def grad to be reached at end of homog inc + crystallite_partionedF0, & ! def grad at start of homog inc + crystallite_subF, & ! def grad to be reached at end of crystallite inc + crystallite_subF0, & ! def grad at start of crystallite inc + crystallite_Lp, & ! current plastic velocitiy grad (end of converged time step) + crystallite_Lp0, & ! plastic velocitiy grad at start of FE inc + crystallite_partionedLp0,& ! plastic velocity grad at start of homog inc + crystallite_subLp0,& ! plastic velocity grad at start of crystallite inc + crystallite_P ! 1st Piola-Kirchhoff stress per grain +real(pReal), dimension (:,:,:,:,:,:,:), allocatable :: crystallite_dPdF, & ! individual dPdF per grain + crystallite_fallbackdPdF ! dPdF fallback for non-converged grains (elastic prediction) -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 +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 CONTAINS @@ -62,7 +61,7 @@ CONTAINS !******************************************************************** ! allocate and initialize per grain variables !******************************************************************** -subroutine crystallite_init() +subroutine crystallite_init(Temperature) !*** variables and functions from other modules ***! use prec, only: pInt, & @@ -83,7 +82,8 @@ subroutine crystallite_init() phase_localConstitution implicit none - !*** input variables ***! + !*** input variables ***! + real(pReal) Temperature !*** output variables ***! @@ -100,13 +100,12 @@ subroutine crystallite_init() iMax = mesh_maxNips eMax = mesh_NcpElems - allocate(crystallite_Temperature(gMax,iMax,eMax)); crystallite_Temperature = 0.0_pReal + 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_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 @@ -139,8 +138,8 @@ 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_partionedTemperature0(g,i,e) = Temperature ! isothermal assumption + 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) @@ -166,7 +165,6 @@ subroutine crystallite_init() 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_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) @@ -281,11 +279,11 @@ subroutine crystallite_stressAndItsTangent(updateJaco) write (6,*) write (6,*) 'Crystallite request from Materialpoint' - 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) + write (6,'(a,/,(f12.7,x))') 'crystallite_partionedTemp0 of 1 1 1' ,crystallite_partionedTemperature0(1,1,1) + write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedF0 of 1 1 1' ,crystallite_partionedF0(1:3,:,1,1,1) + write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedFp0 of 1 1 1' ,crystallite_partionedFp0(1:3,:,1,1,1) + write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedF of 1 1 1' ,crystallite_partionedF(1:3,:,1,1,1) + write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedLp0 of 1 1 1' ,crystallite_partionedLp0(1:3,:,1,1,1) !$OMP PARALLEL DO @@ -456,8 +454,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).AND.& - crystallite_updateTemperature(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 @@ -724,8 +722,7 @@ endsubroutine 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 + crystallite_updateTemperature = abs(residuum/crystallite_Temperature(g,i,e)) < rTol_crystalliteTemperature if (debugger) write(6,'(a,/,f12.4)') 'updated temperature: ', crystallite_Temperature(g,i,e) diff --git a/code/homogenization.f90 b/code/homogenization.f90 index d706ae69b..af939344b 100644 --- a/code/homogenization.f90 +++ b/code/homogenization.f90 @@ -51,7 +51,7 @@ CONTAINS !************************************** !* Module initialization * !************************************** -subroutine homogenization_init() +subroutine homogenization_init(Temperature) use prec, only: pReal,pInt use math, only: math_I3 use IO, only: IO_error, IO_open_file @@ -62,6 +62,7 @@ subroutine homogenization_init() use homogenization_isostrain ! use homogenization_RGC + real(pReal) Temperature integer(pInt), parameter :: fileunit = 200 integer(pInt) e,i,g,myInstance @@ -83,7 +84,7 @@ subroutine homogenization_init() allocate(materialpoint_subF0(3,3,mesh_maxNips,mesh_NcpElems)); materialpoint_subF0 = 0.0_pReal allocate(materialpoint_subF(3,3,mesh_maxNips,mesh_NcpElems)); materialpoint_subF = 0.0_pReal allocate(materialpoint_P(3,3,mesh_maxNips,mesh_NcpElems)); materialpoint_P = 0.0_pReal - allocate(materialpoint_Temperature(mesh_maxNips,mesh_NcpElems)); materialpoint_Temperature = 0.0_pReal + allocate(materialpoint_Temperature(mesh_maxNips,mesh_NcpElems)); materialpoint_Temperature = Temperature allocate(materialpoint_subFrac(mesh_maxNips,mesh_NcpElems)); materialpoint_subFrac = 0.0_pReal allocate(materialpoint_subStep(mesh_maxNips,mesh_NcpElems)); materialpoint_subStep = 0.0_pReal allocate(materialpoint_subdt(mesh_maxNips,mesh_NcpElems)); materialpoint_subdt = 0.0_pReal @@ -180,8 +181,7 @@ subroutine materialpoint_stressAndItsTangent(& use constitutive, only: constitutive_state0, & constitutive_partionedState0, & constitutive_state - use crystallite, only: crystallite_Temperature0, & - crystallite_Temperature, & + use crystallite, only: crystallite_Temperature, & crystallite_F0, & crystallite_Fp0, & crystallite_Fp, & @@ -210,7 +210,8 @@ subroutine materialpoint_stressAndItsTangent(& write (6,*) write (6,*) 'Material Point start' - write (6,'(a,/,3(3(f12.7,x)/))') 'F0 of 8 1',materialpoint_F0(1:3,:,8,1) + write (6,'(a,/,(f12.7,x))') 'Temp0 of 8 1' ,materialpoint_Temperature(8,1) + write (6,'(a,/,3(3(f12.7,x)/))') 'F0 of 8 1',materialpoint_F0(1:3,:,8,1) write (6,'(a,/,3(3(f12.7,x)/))') 'F of 8 1',materialpoint_F(1:3,:,8,1) write (6,'(a,/,3(3(f12.7,x)/))') 'Fp0 of 1 8 1',crystallite_Fp0(1:3,:,1,8,1) write (6,'(a,/,3(3(f12.7,x)/))') 'Lp0 of 1 8 1',crystallite_Lp0(1:3,:,1,8,1) @@ -222,7 +223,7 @@ subroutine materialpoint_stressAndItsTangent(& ! initialize restoration points of grain... 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_partionedTemperature0(1:myNgrains,i,e) = materialpoint_Temperature(i,e) ! ...temperatures crystallite_partionedFp0(:,:,1:myNgrains,i,e) = crystallite_Fp0(:,:,1:myNgrains,i,e) ! ...plastic def grads crystallite_partionedLp0(:,:,1:myNgrains,i,e) = crystallite_Lp0(:,:,1:myNgrains,i,e) ! ...plastic velocity grads crystallite_partionedF0(:,:,1:myNgrains,i,e) = crystallite_F0(:,:,1:myNgrains,i,e) ! ...def grads diff --git a/code/homogenization_isostrain.f90 b/code/homogenization_isostrain.f90 index 64eea767d..009ca6101 100644 --- a/code/homogenization_isostrain.f90 +++ b/code/homogenization_isostrain.f90 @@ -257,7 +257,7 @@ function homogenization_isostrain_averageTemperature(& ! homID = homogenization_typeInstance(mesh_element(3,el)) Ngrains = homogenization_Ngrains(mesh_element(3,el)) - homogenization_isostrain_averageTemperature = sum(Temperature)/Ngrains + homogenization_isostrain_averageTemperature = sum(Temperature(1:Ngrains))/Ngrains return diff --git a/code/numerics.f90 b/code/numerics.f90 index 9754e806c..8457346f2 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -16,7 +16,8 @@ integer(pInt) iJacoStiffness, & ! freque real(pReal) relevantStrain, & ! strain increment considered significant pert_Fg, & ! strain perturbation for FEM Jacobi subStepMin, & ! minimum (relative) size of sub-step allowed during cutback in crystallite - rTol_crystalliteState, & ! relative tolerance in crystallite state loop + rTol_crystalliteState, & ! relative tolerance in crystallite state loop + rTol_crystalliteTemperature, & ! relative tolerance in crystallite temperature loop rTol_crystalliteStress, & ! relative tolerance in crystallite stress loop aTol_crystalliteStress, & ! absolute tolerance in crystallite stress loop resToler, & ! relative tolerance of residual in GIA iteration @@ -88,7 +89,8 @@ subroutine numerics_init() nState = 10_pInt nStress = 40_pInt subStepMin = 1.0e-3_pReal - rTol_crystalliteState = 1.0e-6_pReal + rTol_crystalliteState = 1.0e-6_pReal + rTol_crystalliteTemperature = 1.0e-6_pReal rTol_crystalliteStress = 1.0e-6_pReal aTol_crystalliteStress = 1.0e-8_pReal resToler = 1.0e-4_pReal @@ -129,7 +131,9 @@ subroutine numerics_init() case ('substepmin') subStepMin = IO_floatValue(line,positions,2) case ('rtol_crystallitestate') - rTol_crystalliteState = IO_floatValue(line,positions,2) + rTol_crystalliteState = IO_floatValue(line,positions,2) + case ('rtol_crystalliteTemperature') + rTol_crystalliteTemperature = IO_floatValue(line,positions,2) case ('rtol_crystallitestress') rTol_crystalliteStress = IO_floatValue(line,positions,2) case ('atol_crystallitestress') @@ -164,7 +168,8 @@ subroutine numerics_init() write(6,'(a24,x,i8)') 'nState: ',nState write(6,'(a24,x,i8)') 'nStress: ',nStress write(6,'(a24,x,e8.1)') 'subStepMin: ',subStepMin - write(6,'(a24,x,e8.1)') 'rTol_crystalliteState: ',rTol_crystalliteState + write(6,'(a24,x,e8.1)') 'rTol_crystalliteState: ',rTol_crystalliteState + write(6,'(a24,x,e8.1)') 'rTol_crystalliteTemp: ',rTol_crystalliteTemperature write(6,'(a24,x,e8.1)') 'rTol_crystalliteStress: ',rTol_crystalliteStress write(6,'(a24,x,e8.1)') 'aTol_crystalliteStress: ',aTol_crystalliteStress write(6,'(a24,x,e8.1)') 'resToler: ',resToler @@ -173,7 +178,7 @@ subroutine numerics_init() write(6,'(a24,x,i8)') 'NRiterMax: ',NRiterMax write(6,*) - ! sanity check + ! sanity check (Temperature check missing!!!!!!!) if (relevantStrain <= 0.0_pReal) call IO_error(260) if (iJacoStiffness < 1_pInt) call IO_error(261) if (iJacoLpresiduum < 1_pInt) call IO_error(262)