223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

Major update: corrected treatment of temperature

This commit is contained in:
Luc Hantcherli 2009-07-01 10:29:35 +00:00
parent 2e783df5ed
commit a16b8a619d
7 changed files with 403 additions and 229 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
code/CPFEM.f90
View File

@ -286,7 +286,9 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
! return the local stress and the jacobian from storage ! return the local stress and the jacobian from storage
cauchyStress(1:ngens) = CPFEM_cs(1:ngens,IP,cp_en) cauchyStress(1:ngens) = CPFEM_cs(1:ngens,IP,cp_en)
jacobian(1:ngens,1:ngens) = CPFEM_dcsdE(1:ngens,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 return
end subroutine end subroutine

40
code/constitutive.f90
View File

@ -29,6 +29,7 @@ CONTAINS
!* - constitutive_microstructure !* - constitutive_microstructure
!* - constitutive_LpAndItsTangent !* - constitutive_LpAndItsTangent
!* - constitutive_dotState !* - constitutive_dotState
!* - constitutive_dotTemperature
!* - constitutive_postResults !* - constitutive_postResults
!**************************************** !****************************************
@ -275,6 +276,45 @@ function constitutive_dotState(Tstar_v,Temperature,ipc,ip,el)
end function 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) pure function constitutive_postResults(Tstar_v,Temperature,dt,ipc,ip,el)
!********************************************************************* !*********************************************************************
!* return array of constitutive results * !* return array of constitutive results *

73
code/constitutive_dislobased.f90
View File

@ -93,6 +93,7 @@ CONTAINS
!* - constitutive_microstructure !* - constitutive_microstructure
!* - constitutive_LpAndItsTangent !* - constitutive_LpAndItsTangent
!* - consistutive_dotState !* - consistutive_dotState
!* - constitutive_dotTemperature
!* - consistutive_postResults !* - consistutive_postResults
!**************************************** !****************************************
@ -372,6 +373,7 @@ subroutine constitutive_dislobased_microstructure(Temperature,state,ipc,ip,el)
real(pReal) Temperature real(pReal) Temperature
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state
Temperature = 298.0
matID = phase_constitutionInstance(material_phase(ipc,ip,el)) matID = phase_constitutionInstance(material_phase(ipc,ip,el))
n = constitutive_dislobased_Nslip(matID) n = constitutive_dislobased_Nslip(matID)
!* Quantities derived from state - slip !* 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 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) = & state(ipc,ip,el)%p(6*n+i) = &
(2.0_pReal*kB*Temperature*sqrt(state(ipc,ip,el)%p(2*n+i)))/& (2.0_pReal*kB*Temperature)/(constitutive_dislobased_c1(matID)*constitutive_dislobased_c2(matID)*&
(constitutive_dislobased_c1(matID)*constitutive_dislobased_c3(matID)*constitutive_dislobased_Gmod(matID)*& constitutive_dislobased_c3(matID)*constitutive_dislobased_Gmod(matID)*constitutive_dislobased_bg(matID)**3.0_pReal)*&
state(ipc,ip,el)%p(4*n+i)*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(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)*& 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)) 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) * !* - dLp_dTstar : derivative of Lp (4th-rank tensor) *
!********************************************************************* !*********************************************************************
use prec, only: pReal,pInt,p_vec 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 lattice, only: lattice_Sslip,lattice_Sslip_v
use mesh, only: mesh_NcpElems,mesh_maxNips use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance 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)))) :: & real(pReal), dimension(constitutive_dislobased_Nslip(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: &
gdot_slip,dgdot_dtauslip,tau_slip gdot_slip,dgdot_dtauslip,tau_slip
Temperature = 298.0
matID = phase_constitutionInstance(material_phase(ipc,ip,el)) matID = phase_constitutionInstance(material_phase(ipc,ip,el))
n = constitutive_dislobased_Nslip(matID) n = constitutive_dislobased_Nslip(matID)
@ -453,10 +456,11 @@ subroutine constitutive_dislobased_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Tempera
Lp = 0.0_pReal Lp = 0.0_pReal
gdot_slip = 0.0_pReal gdot_slip = 0.0_pReal
do i = 1,constitutive_dislobased_Nslip(matID) do i = 1,constitutive_dislobased_Nslip(matID)
tau_slip(i) = math_mul6x6(Tstar_v,lattice_Sslip_v(:,i,constitutive_dislobased_structure(matID))) 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) &
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))*& 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)) Lp = Lp + gdot_slip(i)*lattice_Sslip(:,:,i,constitutive_dislobased_structure(matID))
enddo enddo
@ -466,9 +470,11 @@ subroutine constitutive_dislobased_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Tempera
dLp_dTstar = 0.0_pReal dLp_dTstar = 0.0_pReal
dgdot_dtauslip = 0.0_pReal dgdot_dtauslip = 0.0_pReal
do i = 1,constitutive_dislobased_Nslip(matID) do i = 1,constitutive_dislobased_Nslip(matID)
if ((abs(tau_slip(i))-state(ipc,ip,el)%p(3*n+i))>0) & 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)*& 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) & 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) + & dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
dgdot_dtauslip(i)*lattice_Sslip(k,l,i,constitutive_dislobased_structure(matID))* & 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)))) :: & real(pReal), dimension(constitutive_dislobased_Nslip(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: &
constitutive_dislobased_dotState constitutive_dislobased_dotState
Temperature = 298.0
matID = phase_constitutionInstance(material_phase(ipc,ip,el)) matID = phase_constitutionInstance(material_phase(ipc,ip,el))
n = constitutive_dislobased_Nslip(matID) n = constitutive_dislobased_Nslip(matID)
!* Dislocation density evolution !* Dislocation density evolution
constitutive_dislobased_dotState = 0.0_pReal constitutive_dislobased_dotState = 0.0_pReal
do i = 1,n do i = 1,n
tau_slip = dot_product(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)) then 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)*& 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)) 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) 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 constitutive_dislobased_dotState(i) = locks - athermal_recovery
endif endif
enddo enddo
@ -533,6 +538,38 @@ function constitutive_dislobased_dotState(Tstar_v,Temperature,state,ipc,ip,el)
return return
end function 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) pure function constitutive_dislobased_postResults(Tstar_v,Temperature,dt,state,ipc,ip,el)
!********************************************************************* !*********************************************************************
!* return array of constitutive results * !* return array of constitutive results *
@ -544,7 +581,6 @@ pure function constitutive_dislobased_postResults(Tstar_v,Temperature,dt,state,i
!* - el : current element * !* - el : current element *
!********************************************************************* !*********************************************************************
use prec, only: pReal,pInt,p_vec use prec, only: pReal,pInt,p_vec
use math, only: math_mul6x6
use lattice, only: lattice_Sslip_v use lattice, only: lattice_Sslip_v
use mesh, only: mesh_NcpElems,mesh_maxNips use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance,phase_Noutput 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 real(pReal), dimension(6), intent(in) :: Tstar_v
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state
integer(pInt) matID,o,i,c,n 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)))) :: & real(pReal), dimension(constitutive_dislobased_sizePostResults(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: &
constitutive_dislobased_postResults 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)) do o = 1,phase_Noutput(material_phase(ipc,ip,el))
select case(constitutive_dislobased_output(o,matID)) select case(constitutive_dislobased_output(o,matID))
case ('dislodensity') 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(1:n)
c = c + n c = c + n
case ('rateofshear') case ('rateofshear')
do i = 1,n 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 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)*& 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 else
constitutive_dislobased_postResults(c+i) = 0.0_pReal constitutive_dislobased_postResults(c+i) = 0.0_pReal
endif endif
enddo enddo
c = c + n c = c + n
end select end select
enddo enddo

352
code/crystallite.f90
View File

@ -25,7 +25,10 @@ real(pReal), dimension (:,:,:), allocatable :: crystallite_dt, &
crystallite_subdt, & ! substepped time increment of each grain crystallite_subdt, & ! substepped time increment of each grain
crystallite_subFrac, & ! already calculated fraction of increment crystallite_subFrac, & ! already calculated fraction of increment
crystallite_subStep, & ! size of next integration step 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) 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_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_partionedTstar0_v, & ! 2nd Piola-Kirchhoff stress vector at start of homog inc
@ -93,72 +96,39 @@ subroutine crystallite_init()
eMax, & ! maximum number of elements eMax, & ! maximum number of elements
myNgrains 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 gMax = homogenization_maxNgrains
iMax = mesh_maxNips iMax = mesh_maxNips
eMax = mesh_NcpElems 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_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_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_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_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_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_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_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_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_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_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_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_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_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_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_subTstar0_v(6,gMax,iMax,eMax)); crystallite_subTstar0_v = 0.0_pReal
allocate(crystallite_Tstar_v(6,gMax,iMax,eMax)); crystallite_Tstar_v = 0.0_pReal
allocate(crystallite_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_dPdF(3,3,3,3,gMax,iMax,eMax)); crystallite_dPdF = 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_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_dt(gMax,iMax,eMax)); crystallite_dt = 0.0_pReal
allocate(crystallite_subdt(gMax,iMax,eMax)); crystallite_subdt = 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_subFrac(gMax,iMax,eMax)); crystallite_subFrac = 0.0_pReal
allocate(crystallite_subStep(gMax,iMax,eMax)); crystallite_subStep = 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_localConstitution(gMax,iMax,eMax));
allocate(crystallite_requested(gMax,iMax,eMax)); crystallite_requested = .false. allocate(crystallite_requested(gMax,iMax,eMax)); crystallite_requested = .false.
allocate(crystallite_onTrack(gMax,iMax,eMax)); crystallite_onTrack = .false. allocate(crystallite_onTrack(gMax,iMax,eMax)); crystallite_onTrack = .false.
@ -169,14 +139,16 @@ subroutine crystallite_init()
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element
do g = 1,myNgrains 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_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_F0(:,:,g,i,e) = math_I3
crystallite_partionedFp0(:,:,g,i,e) = crystallite_Fp0(:,:,g,i,e) crystallite_partionedFp0(:,:,g,i,e) = crystallite_Fp0(:,:,g,i,e)
crystallite_partionedF0(:,:,g,i,e) = crystallite_F0(:,:,g,i,e) crystallite_partionedF0(:,:,g,i,e) = crystallite_F0(:,:,g,i,e)
crystallite_partionedF(:,:,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_requested(g,i,e) = .true.
crystallite_localConstitution(g,i,e) = phase_localConstitution(material_phase(g,i,e)) crystallite_localConstitution(g,i,e) = phase_localConstitution(material_phase(g,i,e))
enddo
enddo
enddo enddo
enddo enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
@ -190,17 +162,21 @@ subroutine crystallite_init()
write(6,*) '<<<+- crystallite init -+>>>' write(6,*) '<<<+- crystallite init -+>>>'
write(6,*) write(6,*)
write(6,'(a32,x,7(i5,x))') 'crystallite_Nresults: ', crystallite_Nresults 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_Fp: ', shape(crystallite_Fp)
write(6,'(a32,x,7(i5,x))') 'crystallite_Lp: ', shape(crystallite_Lp) 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_Fp0: ', shape(crystallite_Fp0)
write(6,'(a32,x,7(i5,x))') 'crystallite_Lp0: ', shape(crystallite_Lp0) 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_partionedF0: ', shape(crystallite_partionedF0)
write(6,'(a32,x,7(i5,x))') 'crystallite_partionedFp0: ', shape(crystallite_partionedFp0) 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_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_subF0: ', shape(crystallite_subF0)
write(6,'(a32,x,7(i5,x))') 'crystallite_subFp0: ', shape(crystallite_subFp0) 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_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_subdt: ', shape(crystallite_subdt)
write(6,'(a32,x,7(i5,x))') 'crystallite_subFrac: ', shape(crystallite_subFrac) 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_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_requested: ', shape(crystallite_requested)
write(6,'(a32,x,7(i5,x))') 'crystallite_onTrack: ', shape(crystallite_onTrack) 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,'(a32,x,7(i5,x))') 'crystallite_converged: ', shape(crystallite_converged)
@ -277,7 +252,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
!*** output variables ***! !*** output variables ***!
!*** local 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 Fe_guess, & ! guess for elastic deformation gradient
Tstar, & ! 2nd Piola-Kirchhoff stress tensor Tstar, & ! 2nd Piola-Kirchhoff stress tensor
myF, & ! local copy of the deformation gradient 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 logical onTrack, & ! flag indicating wether we are still on track
converged ! flag indicating if iteration converged 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 ------ ! ------ initialize to starting condition ------
write (6,*) write (6,*)
write (6,*) 'Crystallite request from Materialpoint' 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_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_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,/,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 i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains do g = 1,myNgrains
if (crystallite_requested(g,i,e)) then ! initialize restoration point of ... 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 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_subFp0(:,:,g,i,e) = crystallite_partionedFp0(:,:,g,i,e) ! ...plastic def grad
crystallite_subLp0(:,:,g,i,e) = crystallite_partionedLp0(:,:,g,i,e) ! ...plastic velocity grad crystallite_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_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)) 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 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_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 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 constitutive_subState0(g,i,e)%p = constitutive_state(g,i,e)%p ! ...microstructure
@ -405,7 +352,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
endif endif
else else
crystallite_subStep(g,i,e) = 0.5_pReal * crystallite_subStep(g,i,e) ! cut step in half and restore... 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 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 constitutive_state(g,i,e)%p = constitutive_subState0(g,i,e)%p ! ...microstructure
crystallite_Tstar_v(:,g,i,e) = crystallite_subTstar0_v(:,g,i,e) ! ...2nd PK stress crystallite_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) & if ( crystallite_requested(g,i,e) &
.and. crystallite_onTrack(g,i,e) & .and. crystallite_onTrack(g,i,e) &
.and. .not. crystallite_converged(g,i,e)) then ! all undone crystallites .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 crystallite_converged(g,i,e) = .false. ! force at least one iteration step even if state already converged
endif endif
enddo enddo
@ -507,7 +456,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
if ( crystallite_requested(g,i,e) & if ( crystallite_requested(g,i,e) &
.and. crystallite_onTrack(g,i,e) & .and. crystallite_onTrack(g,i,e) &
.and. .not. crystallite_converged(g,i,e)) then ! all undone crystallites .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 if (crystallite_converged(g,i,e)) then
!$OMP CRITICAL (distributionState) !$OMP CRITICAL (distributionState)
debug_StateLoopDistribution(NiterationState) = debug_StateLoopDistribution(NiterationState) + 1 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 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... myState(1:mySizeState) = constitutive_state(g,i,e)%p ! remember unperturbed, converged state...
myF = crystallite_subF(:,:,g,i,e) ! ... and kinematics 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) myFe = crystallite_Fe(:,:,g,i,e)
myLp = crystallite_Lp(:,:,g,i,e) myLp = crystallite_Lp(:,:,g,i,e)
myTstar_v = crystallite_Tstar_v(:,g,i,e) myTstar_v = crystallite_Tstar_v(:,g,i,e)
@ -598,7 +549,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
NiterationState = NiterationState + 1_pInt NiterationState = NiterationState + 1_pInt
if (debugger) write (6,'(a4,x,i6)') 'loop',NiterationState 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) 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 if (debugger) then
write (6,*) '-------------' write (6,*) '-------------'
write (6,'(l,x,l)') onTrack,converged write (6,'(l,x,l)') onTrack,converged
@ -611,7 +563,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
if (converged) & ! converged state warrants stiffness update 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 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... 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_Fe(:,:,g,i,e) = myFe
crystallite_Lp(:,:,g,i,e) = myLp crystallite_Lp(:,:,g,i,e) = myLp
crystallite_Tstar_v(:,g,i,e) = myTstar_v crystallite_Tstar_v(:,g,i,e) = myTstar_v
@ -640,82 +593,147 @@ endsubroutine
!******************************************************************** !********************************************************************
! update the internal state of the constitutive law ! update the internal state of the constitutive law
! and tell whether state has converged ! and tell whether state has converged
!******************************************************************** !********************************************************************
function crystallite_updateState(& function crystallite_updateState(&
g,& ! grain number g,& ! grain number
i,& ! integration point number i,& ! integration point number
e & ! element number e & ! element number
) )
!*** variables and functions from other modules ***! !*** variables and functions from other modules ***!
use prec, only: pReal, & use prec, only: pReal, &
pInt, & pInt, &
pLongInt pLongInt
use numerics, only: rTol_crystalliteState use numerics, only: rTol_crystalliteState
use constitutive, only: constitutive_dotState, & use constitutive, only: constitutive_dotState, &
constitutive_sizeDotState, & constitutive_sizeDotState, &
constitutive_subState0, & constitutive_subState0, &
constitutive_state constitutive_state
use debug, only: debugger, & use debug, only: debugger, &
debug_cumDotStateCalls, & debug_cumDotStateCalls, &
debug_cumDotStateTicks debug_cumDotStateTicks
!*** input variables ***! !*** input variables ***!
integer(pInt), intent(in):: e, & ! element index integer(pInt), intent(in):: e, & ! element index
i, & ! integration point index i, & ! integration point index
g ! grain index g ! grain index
!*** output variables ***! !*** output variables ***!
logical crystallite_updateState ! flag indicating if integration suceeded logical crystallite_updateState ! flag indicating if integration suceeded
!*** local variables ***! !*** local variables ***!
real(pReal), dimension(constitutive_sizeDotState(g,i,e)) :: residuum ! residuum from evolution of microstructure real(pReal), dimension(constitutive_sizeDotState(g,i,e)) :: residuum ! residuum from evolution of microstructure
integer(pInt) mySize integer(pInt) mySize
integer(pLongInt) tick, & integer(pLongInt) tick, &
tock, & tock, &
tickrate, & tickrate, &
maxticks maxticks
!*** global variables ***! !*** global variables ***!
! crystallite_Tstar_v ! crystallite_Tstar_v
! crystallite_subdt ! crystallite_subdt
! crystallite_Temperature ! crystallite_Temperature
mySize = constitutive_sizeDotState(g,i,e) mySize = constitutive_sizeDotState(g,i,e)
! calculate the residuum ! calculate the residuum
call system_clock(count=tick,count_rate=tickrate,count_max=maxticks) 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) - & 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(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) call system_clock(count=tock,count_rate=tickrate,count_max=maxticks)
debug_cumDotStateCalls = debug_cumDotStateCalls + 1_pInt debug_cumDotStateCalls = debug_cumDotStateCalls + 1_pInt
debug_cumDotStateTicks = debug_cumDotStateTicks + tock-tick debug_cumDotStateTicks = debug_cumDotStateTicks + tock-tick
if (tock < tick) debug_cumDotStateTicks = debug_cumDotStateTicks + maxticks if (tock < tick) debug_cumDotStateTicks = debug_cumDotStateTicks + maxticks
! if NaN occured then return without changing the state ! if NaN occured then return without changing the state
if (any(residuum/=residuum)) then if (any(residuum/=residuum)) then
crystallite_updateState = .false. ! indicate state update failed crystallite_updateState = .false. ! indicate state update failed
if (debugger) write(6,*) '::: updateState encountered NaN' if (debugger) write(6,*) '::: updateState encountered NaN'
return return
endif endif
! update the microstructure ! update the microstructure
constitutive_state(g,i,e)%p(1:mySize) = constitutive_state(g,i,e)%p(1:mySize) - residuum 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 ! 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)), & 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 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) if (debugger) write(6,'(a,/,f12.4)') 'updated state: ', constitutive_state(g,i,e)%p(1)
return return
endfunction 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 *** !*** calculation of stress (P) with time integration ***

12
code/debug.f90
View File

@ -11,8 +11,10 @@
integer(pInt), dimension(:), allocatable :: debug_CrystalliteLoopDistribution integer(pInt), dimension(:), allocatable :: debug_CrystalliteLoopDistribution
integer(pLongInt) :: debug_cumLpTicks = 0_pInt integer(pLongInt) :: debug_cumLpTicks = 0_pInt
integer(pLongInt) :: debug_cumDotStateTicks = 0_pInt integer(pLongInt) :: debug_cumDotStateTicks = 0_pInt
integer(pLongInt) :: debug_cumDotTemperatureTicks = 0_pInt
integer(pInt) :: debug_cumLpCalls = 0_pInt integer(pInt) :: debug_cumLpCalls = 0_pInt
integer(pInt) :: debug_cumDotStateCalls = 0_pInt integer(pInt) :: debug_cumDotStateCalls = 0_pInt
integer(pInt) :: debug_cumDotTemperatureCalls = 0_pInt
logical :: debugger = .false. logical :: debugger = .false.
logical :: distribution_init = .false. logical :: distribution_init = .false.
@ -50,8 +52,10 @@ subroutine debug_reset()
debug_CrystalliteLoopDistribution = 0_pInt debug_CrystalliteLoopDistribution = 0_pInt
debug_cumLpTicks = 0_pInt debug_cumLpTicks = 0_pInt
debug_cumDotStateTicks = 0_pInt debug_cumDotStateTicks = 0_pInt
debug_cumDotTemperatureTicks = 0_pInt
debug_cumLpCalls = 0_pInt debug_cumLpCalls = 0_pInt
debug_cumDotStateCalls = 0_pInt debug_cumDotStateCalls = 0_pInt
debug_cumDotTemperatureCalls = 0_pInt
endsubroutine endsubroutine
@ -87,6 +91,14 @@ endsubroutine
dble(debug_cumDotStateTicks)/tickrate/1.0e-6_pReal/debug_cumDotStateCalls dble(debug_cumDotStateTicks)/tickrate/1.0e-6_pReal/debug_cumDotStateCalls
write(6,'(a33,x,i12)') 'total CPU ticks :',debug_cumDotStateTicks write(6,'(a33,x,i12)') 'total CPU ticks :',debug_cumDotStateTicks
endif 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 integral = 0_pInt
write(6,*) write(6,*)

120
code/homogenization.f90
View File

@ -180,13 +180,16 @@ subroutine materialpoint_stressAndItsTangent(&
use constitutive, only: constitutive_state0, & use constitutive, only: constitutive_state0, &
constitutive_partionedState0, & constitutive_partionedState0, &
constitutive_state constitutive_state
use crystallite, only: crystallite_F0, & use crystallite, only: crystallite_Temperature0, &
crystallite_Temperature, &
crystallite_F0, &
crystallite_Fp0, & crystallite_Fp0, &
crystallite_Fp, & crystallite_Fp, &
crystallite_Lp0, & crystallite_Lp0, &
crystallite_Lp, & crystallite_Lp, &
crystallite_Tstar0_v, & crystallite_Tstar0_v, &
crystallite_Tstar_v, & crystallite_Tstar_v, &
crystallite_partionedTemperature0, &
crystallite_partionedF0, & crystallite_partionedF0, &
crystallite_partionedF, & crystallite_partionedF, &
crystallite_partionedFp0, & 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 do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
! initialize restoration points of grain... ! initialize restoration points of grain...
forall (g = 1:myNgrains) constitutive_partionedState0(g,i,e)%p = constitutive_state0(g,i,e)%p ! ...microstructures 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_partionedTemperature0(1:myNgrains,i,e) = crystallite_Temperature0(1:myNgrains,i,e)! ...temperatures
crystallite_partionedLp0(:,:,1:myNgrains,i,e) = crystallite_Lp0(:,:,1:myNgrains,i,e) ! ...plastic velocity grads crystallite_partionedFp0(:,:,1:myNgrains,i,e) = crystallite_Fp0(:,:,1:myNgrains,i,e) ! ...plastic def grads
crystallite_partionedF0(:,:,1:myNgrains,i,e) = crystallite_F0(:,:,1:myNgrains,i,e) ! ...def grads crystallite_partionedLp0(:,:,1:myNgrains,i,e) = crystallite_Lp0(:,:,1:myNgrains,i,e) ! ...plastic velocity grads
crystallite_partionedTstar0_v(:,1:myNgrains,i,e)= crystallite_Tstar0_v(:,1:myNgrains,i,e) ! ...2nd PK stress 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 ... ! initialize restoration points of ...
if (homogenization_sizeState(i,e) > 0_pInt) & if (homogenization_sizeState(i,e) > 0_pInt) &
@ -258,7 +262,8 @@ subroutine materialpoint_stressAndItsTangent(&
if (materialpoint_subStep(i,e) > subStepMin) then if (materialpoint_subStep(i,e) > subStepMin) then
! wind forward grain starting point of... ! 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_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_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 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) materialpoint_subStep(i,e) = 0.5_pReal * materialpoint_subStep(i,e)
! restore... ! 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_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 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 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 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. & if ( materialpoint_requested(i,e) .and. &
.not. materialpoint_doneAndHappy(1,i,e)) then .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 materialpoint_converged(i,e) = all(materialpoint_doneAndHappy(:,i,e)) ! converged if done and happy
endif endif
enddo enddo
@ -358,7 +364,8 @@ subroutine materialpoint_stressAndItsTangent(&
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
call homogenization_averageStressAndItsTangent(i,e) call homogenization_averageStressAndItsTangent(i,e)
call homogenization_averageTemperature(i,e)
enddo enddo
enddo enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
@ -476,38 +483,65 @@ function homogenization_updateState(&
endfunction endfunction
!******************************************************************** !********************************************************************
! derive average stress and stiffness from constituent quantities ! derive average stress and stiffness from constituent quantities
!******************************************************************** !********************************************************************
subroutine homogenization_averageStressAndItsTangent(& subroutine homogenization_averageStressAndItsTangent(&
ip, & ! integration point ip, & ! integration point
el & ! element el & ! element
) )
use prec, only: pReal,pInt use prec, only: pReal,pInt
use mesh, only: mesh_element use mesh, only: mesh_element
use material, only: homogenization_type, homogenization_maxNgrains use material, only: homogenization_type, homogenization_maxNgrains
use crystallite, only: crystallite_P,crystallite_dPdF use crystallite, only: crystallite_P,crystallite_dPdF
use homogenization_isostrain use homogenization_isostrain
implicit none implicit none
integer(pInt), intent(in) :: ip,el integer(pInt), intent(in) :: ip,el
select case(homogenization_type(mesh_element(3,el))) select case(homogenization_type(mesh_element(3,el)))
case (homogenization_isostrain_label) case (homogenization_isostrain_label)
call homogenization_isostrain_averageStressAndItsTangent( materialpoint_P(:,:,ip,el), & call homogenization_isostrain_averageStressAndItsTangent( materialpoint_P(:,:,ip,el), &
materialpoint_dPdF(:,:,:,:,ip,el),& materialpoint_dPdF(:,:,:,:,ip,el),&
crystallite_P(:,:,:,ip,el), & crystallite_P(:,:,:,ip,el), &
crystallite_dPdF(:,:,:,:,:,ip,el), & crystallite_dPdF(:,:,:,:,:,ip,el), &
ip, & ip, &
el) el)
end select end select
return return
endsubroutine 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 ! return array of homogenization results for post file inclusion
! call only, if homogenization_sizePostResults(ip,el) > 0 !! ! call only, if homogenization_sizePostResults(ip,el) > 0 !!

31
code/homogenization_isostrain.f90
View File

@ -232,7 +232,36 @@ subroutine homogenization_isostrain_averageStressAndItsTangent(&
return 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
!******************************************************************** !********************************************************************