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

now with first draft of nonlocal constitutive law 

debugging memory leak closed
debugging counters corrected

center of gravity stored in mesh

state updated is now split into a collecting loop and an execution

updateState and updateTemperature fill sequentially separate logicals and evaluate afterwards to converged

added 3x3 transposition function, norm for 3x1 matrix and 33x3 matrix multiplication in math

non-converged crystallite triggers materialpoint cutback (used to respond elastically)

non-converged materialpoint raises terminal illness which in turn renders whole FE increment useless by means of odd stress/stiffness and thus waits for FE cutback
This commit is contained in:
Christoph Kords 2009-08-11 16:31:57 +00:00
parent 4689966c79
commit 8ed3ddc03b
15 changed files with 1695 additions and 359 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

14
code/CPFEM.f90
View File

@ -73,6 +73,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
use FEsolving, only: FE_init, &
parallelExecution, &
outdatedFFN1, &
terminallyIll, &
cycleCounter, &
theInc, &
theCycle, &
@ -201,7 +202,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
j = 1:mesh_maxNips, &
k = 1:mesh_NcpElems ) &
constitutive_state0(i,j,k)%p = constitutive_state(i,j,k)%p ! microstructure of crystallites
write(6,'(a10,/,4(3(f20.8,x),/))') 'aged state',constitutive_state(1,1,1)%p
write(6,'(a10,/,4(3(e20.8,x),/))') 'aged state',constitutive_state(1,1,1)%p
do k = 1,mesh_NcpElems
do j = 1,mesh_maxNips
if (homogenization_sizeState(j,k) > 0_pInt) &
@ -211,10 +212,11 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
endif
! deformation gradient outdated or any actual deformation gradient differs more than relevantStrain from the stored one
if (outdatedFFN1 .or. any(abs(ffn1 - materialpoint_F(:,:,IP,cp_en)) > relevantStrain)) then
if (.not. outdatedFFN1) &
if (terminallyIll .or. outdatedFFN1 .or. any(abs(ffn1 - materialpoint_F(:,:,IP,cp_en)) > relevantStrain)) then
if (.not. terminallyIll .and. .not. outdatedFFN1) then
write(6,'(a11,x,i5,x,i2,x,a10,/,3(3(f10.3,x),/))') 'outdated at',cp_en,IP,'FFN1 now:',ffn1(:,1),ffn1(:,2),ffn1(:,3)
outdatedFFN1 = .true.
endif
CPFEM_cs(1:ngens,IP,cp_en) = CPFEM_odd_stress
CPFEM_dcsde(1:ngens,1:ngens,IP,cp_en) = CPFEM_odd_jacobian*math_identity2nd(ngens)
@ -240,6 +242,10 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
CPFEM_calc_done = .true.
endif
if (terminallyIll) then
CPFEM_cs(1:ngens,IP,cp_en) = CPFEM_odd_stress
CPFEM_dcsde(1:ngens,1:ngens,IP,cp_en) = CPFEM_odd_jacobian*math_identity2nd(ngens)
else
! translate from P to CS
Kirchhoff = math_mul33x33(materialpoint_P(:,:,IP, cp_en),transpose(materialpoint_F(:,:,IP, cp_en)))
J_inverse = 1.0_pReal/math_det3x3(materialpoint_F(:,:,IP, cp_en))
@ -258,7 +264,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
forall(i=1:3,j=1:3,k=1:3,l=1:3) &
H_sym(i,j,k,l)= 0.25_pReal*(H(i,j,k,l)+H(j,i,k,l)+H(i,j,l,k)+H(j,i,l,k)) ! where to use the symmetric version??
CPFEM_dcsde(1:ngens,1:ngens,IP,cp_en) = math_Mandel3333to66(J_inverse*H)
endif
endif
! --+>> COLLECTION OF FEM DATA AND RETURN OF ODD STRESS AND JACOBIAN <<+--

2
code/FEsolving.f90
View File

@ -9,7 +9,7 @@
integer(pInt) cycleCounter
integer(pInt) theInc,theCycle,theLovl
real(pReal) theTime
logical :: lastIncConverged = .false.,outdatedByNewInc = .false.,outdatedFFN1 = .false.
logical :: lastIncConverged = .false.,outdatedByNewInc = .false.,outdatedFFN1 = .false.,terminallyIll = .false.
logical :: symmetricSolver = .false.
logical :: parallelExecution = .true.
integer(pInt), dimension(:,:), allocatable :: FEsolving_execIP

5
code/IO.f90
View File

@ -858,7 +858,7 @@ endfunction
case (265)
msg = 'Limit for crystallite loop too small'
case (266)
msg = 'Limit for state loop too small'
msg = 'Limit for crystallite state loop too small'
case (267)
msg = 'Limit for stress loop too small'
case (268)
@ -884,6 +884,9 @@ endfunction
case (277)
msg = 'Non-positive relevant mismatch in RGC'
case (279)
msg = 'Limit for materialpoint state loop too small'
case (300)
msg = 'This material can only be used with elements with three direct stress components'
case (500)

131
code/constitutive.f90
View File

@ -16,11 +16,14 @@ MODULE constitutive
type(p_vec), dimension(:,:,:), allocatable :: constitutive_state0, & ! pointer array to microstructure at start of FE inc
constitutive_partionedState0, & ! pointer array to microstructure at start of homogenization inc
constitutive_subState0, & ! pointer array to microstructure at start of crystallite inc
constitutive_state ! pointer array to current microstructure (end of converged time step)
constitutive_state, & ! pointer array to current microstructure (end of converged time step)
constitutive_dotState ! pointer array to evolution of current microstructure
integer(pInt), dimension(:,:,:), allocatable :: constitutive_sizeDotState, & ! size of dotState array
constitutive_sizeState, & ! size of state array per grain
constitutive_sizePostResults ! size of postResults array per grain
integer(pInt) constitutive_maxSizeDotState,constitutive_maxSizeState,constitutive_maxSizePostResults
integer(pInt) constitutive_maxSizeDotState, &
constitutive_maxSizeState, &
constitutive_maxSizePostResults
CONTAINS
!****************************************
@ -28,8 +31,8 @@ CONTAINS
!* - constitutive_homogenizedC
!* - constitutive_microstructure
!* - constitutive_LpAndItsTangent
!* - constitutive_dotState
!* - constitutive_dotTemperature
!* - constitutive_collectDotState
!* - constitutive_collectDotTemperature
!* - constitutive_postResults
!****************************************
@ -46,6 +49,7 @@ subroutine constitutive_init()
use constitutive_j2
use constitutive_phenopowerlaw
use constitutive_dislobased
use constitutive_nonlocal
integer(pInt), parameter :: fileunit = 200
integer(pInt) e,i,g,p,myInstance,myNgrains
@ -58,6 +62,7 @@ subroutine constitutive_init()
call constitutive_j2_init(fileunit) ! parse all phases of this constitution
call constitutive_phenopowerlaw_init(fileunit)
call constitutive_dislobased_init(fileunit)
call constitutive_nonlocal_init(fileunit)
close(fileunit)
@ -78,6 +83,9 @@ subroutine constitutive_init()
case (constitutive_dislobased_label)
thisOutput => constitutive_dislobased_output
thisSize => constitutive_dislobased_sizePostResult
case (constitutive_nonlocal_label)
thisOutput => constitutive_nonlocal_output
thisSize => constitutive_nonlocal_sizePostResult
case default
knownConstitution = .false.
end select
@ -101,10 +109,10 @@ subroutine constitutive_init()
allocate(constitutive_partionedState0(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
allocate(constitutive_subState0(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
allocate(constitutive_state(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
allocate(constitutive_dotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
allocate(constitutive_sizeDotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_sizeDotState = 0_pInt
allocate(constitutive_sizeState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_sizeState = 0_pInt
allocate(constitutive_sizePostResults(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_sizePostResults = 0_pInt
allocate(constitutive_sizePostResults(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)); constitutive_sizePostResults = 0_pInt
do e = 1,mesh_NcpElems ! loop over elements
myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = 1,FE_Nips(mesh_element(2,e)) ! loop over IPs
@ -117,28 +125,41 @@ subroutine constitutive_init()
allocate(constitutive_partionedState0(g,i,e)%p(constitutive_j2_sizeState(myInstance)))
allocate(constitutive_subState0(g,i,e)%p(constitutive_j2_sizeState(myInstance)))
allocate(constitutive_state(g,i,e)%p(constitutive_j2_sizeState(myInstance)))
allocate(constitutive_dotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance)))
constitutive_state0(g,i,e)%p = constitutive_j2_stateInit(myInstance)
constitutive_sizeDotState(g,i,e) = constitutive_j2_sizeDotState(myInstance)
constitutive_sizeState(g,i,e) = constitutive_j2_sizeState(myInstance)
constitutive_sizeDotState(g,i,e) = constitutive_j2_sizeDotState(myInstance)
constitutive_sizePostResults(g,i,e) = constitutive_j2_sizePostResults(myInstance)
case (constitutive_phenopowerlaw_label)
allocate(constitutive_state0(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance)))
allocate(constitutive_partionedState0(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance)))
allocate(constitutive_subState0(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance)))
allocate(constitutive_state(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance)))
allocate(constitutive_dotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance)))
constitutive_state0(g,i,e)%p = constitutive_phenopowerlaw_stateInit(myInstance)
constitutive_sizeDotState(g,i,e) = constitutive_phenopowerlaw_sizeDotState(myInstance)
constitutive_sizeState(g,i,e) = constitutive_phenopowerlaw_sizeState(myInstance)
constitutive_sizeDotState(g,i,e) = constitutive_phenopowerlaw_sizeDotState(myInstance)
constitutive_sizePostResults(g,i,e) = constitutive_phenopowerlaw_sizePostResults(myInstance)
case (constitutive_dislobased_label)
allocate(constitutive_state0(g,i,e)%p(constitutive_dislobased_sizeState(myInstance)))
allocate(constitutive_partionedState0(g,i,e)%p(constitutive_dislobased_sizeState(myInstance)))
allocate(constitutive_subState0(g,i,e)%p(constitutive_dislobased_sizeState(myInstance)))
allocate(constitutive_state(g,i,e)%p(constitutive_dislobased_sizeState(myInstance)))
allocate(constitutive_dotState(g,i,e)%p(constitutive_dislobased_sizeDotState(myInstance)))
constitutive_state0(g,i,e)%p = constitutive_dislobased_stateInit(myInstance)
constitutive_sizeDotState(g,i,e) = constitutive_dislobased_sizeDotState(myInstance)
constitutive_sizeState(g,i,e) = constitutive_dislobased_sizeState(myInstance)
constitutive_sizeDotState(g,i,e) = constitutive_dislobased_sizeDotState(myInstance)
constitutive_sizePostResults(g,i,e) = constitutive_dislobased_sizePostResults(myInstance)
case (constitutive_nonlocal_label)
allocate(constitutive_state0(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance)))
allocate(constitutive_partionedState0(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance)))
allocate(constitutive_subState0(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance)))
allocate(constitutive_state(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance)))
allocate(constitutive_dotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance)))
constitutive_state0(g,i,e)%p = constitutive_nonlocal_stateInit(myInstance)
constitutive_sizeState(g,i,e) = constitutive_nonlocal_sizeState(myInstance)
constitutive_sizeDotState(g,i,e) = constitutive_nonlocal_sizeDotState(myInstance)
constitutive_sizePostResults(g,i,e) = constitutive_nonlocal_sizePostResults(myInstance)
case default
call IO_error(200,material_phase(g,i,e)) ! unknown constitution
end select
@ -146,8 +167,9 @@ subroutine constitutive_init()
enddo
enddo
enddo
constitutive_maxSizeDotState = maxval(constitutive_sizeDotState)
constitutive_maxSizeState = maxval(constitutive_sizeState)
constitutive_maxSizeDotState = maxval(constitutive_sizeDotState)
constitutive_maxSizePostResults = maxval(constitutive_sizePostResults)
write(6,*)
@ -157,6 +179,7 @@ subroutine constitutive_init()
write(6,'(a32,x,7(i5,x))') 'constitutive_partionedState0: ', shape(constitutive_partionedState0)
write(6,'(a32,x,7(i5,x))') 'constitutive_subState0: ', shape(constitutive_subState0)
write(6,'(a32,x,7(i5,x))') 'constitutive_state: ', shape(constitutive_state)
write(6,'(a32,x,7(i5,x))') 'constitutive_dotState: ', shape(constitutive_dotState)
write(6,'(a32,x,7(i5,x))') 'constitutive_sizeState: ', shape(constitutive_sizeState)
write(6,'(a32,x,7(i5,x))') 'constitutive_sizeDotState: ', shape(constitutive_sizeDotState)
write(6,'(a32,x,7(i5,x))') 'constitutive_sizePostResults: ', shape(constitutive_sizePostResults)
@ -166,7 +189,7 @@ subroutine constitutive_init()
return
end subroutine
endsubroutine
function constitutive_homogenizedC(ipc,ip,el)
@ -183,6 +206,7 @@ function constitutive_homogenizedC(ipc,ip,el)
use constitutive_j2
use constitutive_phenopowerlaw
use constitutive_dislobased
use constitutive_nonlocal
implicit none
!* Definition of variables
@ -190,19 +214,26 @@ function constitutive_homogenizedC(ipc,ip,el)
real(pReal), dimension(6,6) :: constitutive_homogenizedC
select case (phase_constitution(material_phase(ipc,ip,el)))
case (constitutive_j2_label)
constitutive_homogenizedC = constitutive_j2_homogenizedC(constitutive_state,ipc,ip,el)
case (constitutive_phenopowerlaw_label)
constitutive_homogenizedC = constitutive_phenopowerlaw_homogenizedC(constitutive_state,ipc,ip,el)
case (constitutive_dislobased_label)
constitutive_homogenizedC = constitutive_dislobased_homogenizedC(constitutive_state,ipc,ip,el)
case (constitutive_nonlocal_label)
constitutive_homogenizedC = constitutive_nonlocal_homogenizedC(constitutive_state,ipc,ip,el)
end select
return
end function
endfunction
subroutine constitutive_microstructure(Temperature,ipc,ip,el)
subroutine constitutive_microstructure(Temperature,Fp,ipc,ip,el)
!*********************************************************************
!* This function calculates from state needed variables *
!* INPUT: *
@ -213,29 +244,42 @@ subroutine constitutive_microstructure(Temperature,ipc,ip,el)
!* - el : current element *
!*********************************************************************
use prec, only: pReal,pInt
use material, only: phase_constitution,material_phase
use material, only: phase_constitution, &
material_phase, &
homogenization_maxNgrains
use mesh, only: mesh_NcpElems, &
mesh_maxNips
use constitutive_j2
use constitutive_phenopowerlaw
use constitutive_dislobased
use constitutive_nonlocal
implicit none
!* Definition of variables
integer(pInt) ipc,ip,el
real(pReal) Temperature
integer(pInt), intent(in) :: ipc,ip,el
real(pReal), intent(in) :: Temperature
real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: Fp
select case (phase_constitution(material_phase(ipc,ip,el)))
case (constitutive_j2_label)
call constitutive_j2_microstructure(Temperature,constitutive_state,ipc,ip,el)
case (constitutive_phenopowerlaw_label)
call constitutive_phenopowerlaw_microstructure(Temperature,constitutive_state,ipc,ip,el)
case (constitutive_dislobased_label)
call constitutive_dislobased_microstructure(Temperature,constitutive_state,ipc,ip,el)
case (constitutive_nonlocal_label)
call constitutive_nonlocal_microstructure(Temperature, Fp, constitutive_state,ipc,ip,el)
end select
end subroutine
endsubroutine
subroutine constitutive_LpAndItsTangent(Lp,dLp_dTstar, Tstar_v,Temperature,ipc,ip,el)
subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, Temperature, ipc, ip, el)
!*********************************************************************
!* This subroutine contains the constitutive equation for *
!* calculating the velocity gradient *
@ -253,6 +297,7 @@ subroutine constitutive_LpAndItsTangent(Lp,dLp_dTstar, Tstar_v,Temperature,ipc,i
use constitutive_j2
use constitutive_phenopowerlaw
use constitutive_dislobased
use constitutive_nonlocal
implicit none
!* Definition of variables
@ -263,19 +308,26 @@ subroutine constitutive_LpAndItsTangent(Lp,dLp_dTstar, Tstar_v,Temperature,ipc,i
real(pReal), dimension(9,9) :: dLp_dTstar
select case (phase_constitution(material_phase(ipc,ip,el)))
case (constitutive_j2_label)
call constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_phenopowerlaw_label)
call constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_dislobased_label)
call constitutive_dislobased_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_nonlocal_label)
call constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, Temperature, constitutive_state, ipc, ip, el)
end select
return
end subroutine
endsubroutine
function constitutive_dotState(Tstar_v,Temperature,ipc,ip,el)
subroutine constitutive_collectDotState(Tstar_v, Fp, invFp, Temperature, ipc, ip, el)
!*********************************************************************
!* This subroutine contains the constitutive equation for *
!* calculating the rate of change of microstructure *
@ -289,28 +341,37 @@ function constitutive_dotState(Tstar_v,Temperature,ipc,ip,el)
!* - constitutive_dotState : evolution of state variable *
!*********************************************************************
use prec, only: pReal,pInt
use debug
use material, only: phase_constitution,material_phase
use constitutive_j2
use constitutive_phenopowerlaw
use constitutive_dislobased
use constitutive_nonlocal
implicit none
!* Definition of variables
integer(pInt) ipc,ip,el
real(pReal) Temperature
real(pReal), dimension(3,3) :: Fp, invFp
real(pReal), dimension(6) :: Tstar_v
real(pReal), dimension(constitutive_sizeDotState(ipc,ip,el)) :: constitutive_dotState
select case (phase_constitution(material_phase(ipc,ip,el)))
case (constitutive_j2_label)
constitutive_dotState = constitutive_j2_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el)
constitutive_dotState(ipc,ip,el)%p = constitutive_j2_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_phenopowerlaw_label)
constitutive_dotState = constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el)
constitutive_dotState(ipc,ip,el)%p = constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_dislobased_label)
constitutive_dotState = constitutive_dislobased_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el)
constitutive_dotState(ipc,ip,el)%p = constitutive_dislobased_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_nonlocal_label)
call constitutive_nonlocal_dotState(constitutive_dotState,Tstar_v,Fp,invFp,Temperature,constitutive_state,ipc,ip,el)
end select
return
end function
endsubroutine
function constitutive_dotTemperature(Tstar_v,Temperature,ipc,ip,el)
@ -331,25 +392,32 @@ function constitutive_dotTemperature(Tstar_v,Temperature,ipc,ip,el)
use constitutive_j2
use constitutive_phenopowerlaw
use constitutive_dislobased
use constitutive_nonlocal
implicit none
!* Definition of variables
integer(pInt) ipc,ip,el
real(pReal) Temperature
real(pReal), dimension(6) :: Tstar_v
real(pReal) constitutive_dotTemperature
real(pReal), dimension(6) :: Tstar_v
select case (phase_constitution(material_phase(ipc,ip,el)))
case (constitutive_j2_label)
constitutive_dotTemperature = constitutive_j2_dotTemperature(Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_phenopowerlaw_label)
constitutive_dotTemperature = constitutive_phenopowerlaw_dotTemperature(Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_dislobased_label)
constitutive_dotTemperature = constitutive_dislobased_dotTemperature(Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_nonlocal_label)
constitutive_dotTemperature = constitutive_nonlocal_dotTemperature(Tstar_v,Temperature,constitutive_state,ipc,ip,el)
end select
return
end function
endfunction
pure function constitutive_postResults(Tstar_v,Temperature,dt,ipc,ip,el)
@ -367,6 +435,7 @@ pure function constitutive_postResults(Tstar_v,Temperature,dt,ipc,ip,el)
use constitutive_j2
use constitutive_phenopowerlaw
use constitutive_dislobased
use constitutive_nonlocal
implicit none
!* Definition of variables
@ -377,17 +446,23 @@ pure function constitutive_postResults(Tstar_v,Temperature,dt,ipc,ip,el)
constitutive_postResults = 0.0_pReal
select case (phase_constitution(material_phase(ipc,ip,el)))
case (constitutive_j2_label)
constitutive_postResults = constitutive_j2_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el)
case (constitutive_phenopowerlaw_label)
constitutive_postResults = constitutive_phenopowerlaw_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el)
case (constitutive_dislobased_label)
constitutive_postResults = constitutive_dislobased_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el)
case (constitutive_nonlocal_label)
constitutive_postResults = constitutive_nonlocal_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el)
end select
return
end function
endfunction
END MODULE

1113
code/constitutive_nonlocal.f90 Normal file
View File

File diff suppressed because it is too large Load Diff

231
code/crystallite.f90
View File

@ -34,6 +34,7 @@ real(pReal), dimension (:,:,:,:), allocatable :: crystallite_Tstar_v, &
crystallite_subTstar0_v ! 2nd Piola-Kirchhoff stress vector at start of crystallite inc
real(pReal), dimension (:,:,:,:,:), allocatable :: crystallite_Fe, & ! current "elastic" def grad (end of converged time step)
crystallite_Fp, & ! current plastic def grad (end of converged time step)
crystallite_invFp, & ! inverse of current plastic def grad (end of converged time step)
crystallite_Fp0, & ! plastic def grad at start of FE inc
crystallite_partionedFp0,& ! plastic def grad at start of homog inc
crystallite_subFp0,& ! plastic def grad at start of crystallite inc
@ -53,7 +54,10 @@ real(pReal), dimension (:,:,:,:,:,:,:), allocatable :: crystallite_dPdF, &
logical, dimension (:,:,:), allocatable :: crystallite_localConstitution, & ! indicates this grain to have purely local constitutive law
crystallite_requested, & ! flag to request crystallite calculation
crystallite_onTrack, & ! flag to indicate ongoing calculation
crystallite_converged ! convergence flag
crystallite_converged, & ! convergence flag
crystallite_stateConverged, & ! flag indicating convergence of state
crystallite_temperatureConverged, & ! flag indicating convergence of temperature
crystallite_nonfinished ! requested and ontrack but not converged
CONTAINS
@ -103,6 +107,7 @@ subroutine crystallite_init(Temperature)
allocate(crystallite_P(3,3,gMax,iMax,eMax)); crystallite_P = 0.0_pReal
allocate(crystallite_Fe(3,3,gMax,iMax,eMax)); crystallite_Fe = 0.0_pReal
allocate(crystallite_Fp(3,3,gMax,iMax,eMax)); crystallite_Fp = 0.0_pReal
allocate(crystallite_invFp(3,3,gMax,iMax,eMax)); crystallite_invFp = 0.0_pReal
allocate(crystallite_Lp(3,3,gMax,iMax,eMax)); crystallite_Lp = 0.0_pReal
allocate(crystallite_Tstar_v(6,gMax,iMax,eMax)); crystallite_Tstar_v = 0.0_pReal
allocate(crystallite_F0(3,3,gMax,iMax,eMax)); crystallite_F0 = 0.0_pReal
@ -129,8 +134,11 @@ subroutine crystallite_init(Temperature)
allocate(crystallite_subStep(gMax,iMax,eMax)); crystallite_subStep = 0.0_pReal
allocate(crystallite_localConstitution(gMax,iMax,eMax)); crystallite_localConstitution = .true.
allocate(crystallite_requested(gMax,iMax,eMax)); crystallite_requested = .false.
allocate(crystallite_onTrack(gMax,iMax,eMax)); crystallite_onTrack = .false.
allocate(crystallite_onTrack(gMax,iMax,eMax)); crystallite_onTrack = .true.
allocate(crystallite_converged(gMax,iMax,eMax)); crystallite_converged = .true.
allocate(crystallite_stateConverged(gMax,iMax,eMax)); crystallite_stateConverged = .false.
allocate(crystallite_temperatureConverged(gMax,iMax,eMax)); crystallite_temperatureConverged = .false.
allocate(crystallite_nonfinished(gMax,iMax,eMax)); crystallite_nonfinished = .true.
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over all cp elements
@ -158,39 +166,43 @@ subroutine crystallite_init(Temperature)
write(6,*)
write(6,*) '<<<+- crystallite init -+>>>'
write(6,*)
write(6,'(a32,x,7(i5,x))') 'crystallite_Nresults: ', crystallite_Nresults
write(6,'(a32,x,7(i5,x))') 'crystallite_Temperature ', shape(crystallite_Temperature)
write(6,'(a32,x,7(i5,x))') 'crystallite_Fe: ', shape(crystallite_Fe)
write(6,'(a32,x,7(i5,x))') 'crystallite_Fp: ', shape(crystallite_Fp)
write(6,'(a32,x,7(i5,x))') 'crystallite_Lp: ', shape(crystallite_Lp)
write(6,'(a32,x,7(i5,x))') 'crystallite_F0: ', shape(crystallite_F0)
write(6,'(a32,x,7(i5,x))') 'crystallite_Fp0: ', shape(crystallite_Fp0)
write(6,'(a32,x,7(i5,x))') 'crystallite_Lp0: ', shape(crystallite_Lp0)
write(6,'(a32,x,7(i5,x))') 'crystallite_partionedF: ', shape(crystallite_partionedF)
write(6,'(a32,x,7(i5,x))') 'crystallite_partionedTemp0: ', shape(crystallite_partionedTemperature0)
write(6,'(a32,x,7(i5,x))') 'crystallite_partionedF0: ', shape(crystallite_partionedF0)
write(6,'(a32,x,7(i5,x))') 'crystallite_partionedFp0: ', shape(crystallite_partionedFp0)
write(6,'(a32,x,7(i5,x))') 'crystallite_partionedLp0: ', shape(crystallite_partionedLp0)
write(6,'(a32,x,7(i5,x))') 'crystallite_subF: ', shape(crystallite_subF)
write(6,'(a32,x,7(i5,x))') 'crystallite_subTemperature0: ', shape(crystallite_subTemperature0)
write(6,'(a32,x,7(i5,x))') 'crystallite_subF0: ', shape(crystallite_subF0)
write(6,'(a32,x,7(i5,x))') 'crystallite_subFp0: ', shape(crystallite_subFp0)
write(6,'(a32,x,7(i5,x))') 'crystallite_subLp0: ', shape(crystallite_subLp0)
write(6,'(a32,x,7(i5,x))') 'crystallite_P: ', shape(crystallite_P)
write(6,'(a32,x,7(i5,x))') 'crystallite_Tstar_v: ', shape(crystallite_Tstar_v)
write(6,'(a32,x,7(i5,x))') 'crystallite_Tstar0_v: ', shape(crystallite_Tstar0_v)
write(6,'(a32,x,7(i5,x))') 'crystallite_partionedTstar0_v: ', shape(crystallite_partionedTstar0_v)
write(6,'(a32,x,7(i5,x))') 'crystallite_subTstar0_v: ', shape(crystallite_subTstar0_v)
write(6,'(a32,x,7(i5,x))') 'crystallite_dPdF: ', shape(crystallite_dPdF)
write(6,'(a32,x,7(i5,x))') 'crystallite_fallbackdPdF: ', shape(crystallite_fallbackdPdF)
write(6,'(a32,x,7(i5,x))') 'crystallite_dt: ', shape(crystallite_dt)
write(6,'(a32,x,7(i5,x))') 'crystallite_subdt: ', shape(crystallite_subdt)
write(6,'(a32,x,7(i5,x))') 'crystallite_subFrac: ', shape(crystallite_subFrac)
write(6,'(a32,x,7(i5,x))') 'crystallite_subStep: ', shape(crystallite_subStep)
write(6,'(a32,x,7(i5,x))') 'crystallite_localConstitution: ', shape(crystallite_localConstitution)
write(6,'(a32,x,7(i5,x))') 'crystallite_requested: ', shape(crystallite_requested)
write(6,'(a32,x,7(i5,x))') 'crystallite_onTrack: ', shape(crystallite_onTrack)
write(6,'(a32,x,7(i5,x))') 'crystallite_converged: ', shape(crystallite_converged)
write(6,'(a35,x,7(i5,x))') 'crystallite_Nresults: ', crystallite_Nresults
write(6,*)
write(6,'(a35,x,7(i5,x))') 'crystallite_Temperature: ', shape(crystallite_Temperature)
write(6,'(a35,x,7(i5,x))') 'crystallite_Fe: ', shape(crystallite_Fe)
write(6,'(a35,x,7(i5,x))') 'crystallite_Fp: ', shape(crystallite_Fp)
write(6,'(a35,x,7(i5,x))') 'crystallite_Lp: ', shape(crystallite_Lp)
write(6,'(a35,x,7(i5,x))') 'crystallite_F0: ', shape(crystallite_F0)
write(6,'(a35,x,7(i5,x))') 'crystallite_Fp0: ', shape(crystallite_Fp0)
write(6,'(a35,x,7(i5,x))') 'crystallite_Lp0: ', shape(crystallite_Lp0)
write(6,'(a35,x,7(i5,x))') 'crystallite_partionedF: ', shape(crystallite_partionedF)
write(6,'(a35,x,7(i5,x))') 'crystallite_partionedTemp0: ', shape(crystallite_partionedTemperature0)
write(6,'(a35,x,7(i5,x))') 'crystallite_partionedF0: ', shape(crystallite_partionedF0)
write(6,'(a35,x,7(i5,x))') 'crystallite_partionedFp0: ', shape(crystallite_partionedFp0)
write(6,'(a35,x,7(i5,x))') 'crystallite_partionedLp0: ', shape(crystallite_partionedLp0)
write(6,'(a35,x,7(i5,x))') 'crystallite_subF: ', shape(crystallite_subF)
write(6,'(a35,x,7(i5,x))') 'crystallite_subTemperature0: ', shape(crystallite_subTemperature0)
write(6,'(a35,x,7(i5,x))') 'crystallite_subF0: ', shape(crystallite_subF0)
write(6,'(a35,x,7(i5,x))') 'crystallite_subFp0: ', shape(crystallite_subFp0)
write(6,'(a35,x,7(i5,x))') 'crystallite_subLp0: ', shape(crystallite_subLp0)
write(6,'(a35,x,7(i5,x))') 'crystallite_P: ', shape(crystallite_P)
write(6,'(a35,x,7(i5,x))') 'crystallite_Tstar_v: ', shape(crystallite_Tstar_v)
write(6,'(a35,x,7(i5,x))') 'crystallite_Tstar0_v: ', shape(crystallite_Tstar0_v)
write(6,'(a35,x,7(i5,x))') 'crystallite_partionedTstar0_v: ', shape(crystallite_partionedTstar0_v)
write(6,'(a35,x,7(i5,x))') 'crystallite_subTstar0_v: ', shape(crystallite_subTstar0_v)
write(6,'(a35,x,7(i5,x))') 'crystallite_dPdF: ', shape(crystallite_dPdF)
write(6,'(a35,x,7(i5,x))') 'crystallite_fallbackdPdF: ', shape(crystallite_fallbackdPdF)
write(6,'(a35,x,7(i5,x))') 'crystallite_dt: ', shape(crystallite_dt)
write(6,'(a35,x,7(i5,x))') 'crystallite_subdt: ', shape(crystallite_subdt)
write(6,'(a35,x,7(i5,x))') 'crystallite_subFrac: ', shape(crystallite_subFrac)
write(6,'(a35,x,7(i5,x))') 'crystallite_subStep: ', shape(crystallite_subStep)
write(6,'(a35,x,7(i5,x))') 'crystallite_localConstitution: ', shape(crystallite_localConstitution)
write(6,'(a35,x,7(i5,x))') 'crystallite_requested: ', shape(crystallite_requested)
write(6,'(a35,x,7(i5,x))') 'crystallite_onTrack: ', shape(crystallite_onTrack)
write(6,'(a35,x,7(i5,x))') 'crystallite_converged: ', shape(crystallite_converged)
write(6,'(a35,x,7(i5,x))') 'crystallite_stateConverged: ', shape(crystallite_stateConverged)
write(6,'(a35,x,7(i5,x))') 'crystallite_temperatureConverged: ', shape(crystallite_temperatureConverged)
write(6,'(a35,x,7(i5,x))') 'crystallite_nonfinished: ', shape(crystallite_nonfinished)
write(6,*)
write(6,*) 'Number of non-local grains: ',count(.not. crystallite_localConstitution)
call flush(6)
@ -219,7 +231,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
nCryst
use debug, only: debugger, &
debug_CrystalliteLoopDistribution, &
debug_StateLoopDistribution, &
debug_CrystalliteStateLoopDistribution, &
debug_StiffnessStateLoopDistribution
use IO, only: IO_warning
use math, only: math_inv3x3, &
@ -235,10 +247,14 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
use material, only: homogenization_Ngrains
use constitutive, only: constitutive_maxSizeState, &
constitutive_sizeState, &
constitutive_sizeDotState, &
constitutive_state, &
constitutive_subState0, &
constitutive_partionedState0, &
constitutive_homogenizedC
constitutive_homogenizedC, &
constitutive_dotState, &
constitutive_collectDotState, &
constitutive_dotTemperature
implicit none
@ -258,7 +274,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
myLp, & ! local copy of the plastic velocity gradient
myP ! local copy of the 1st Piola-Kirchhoff stress tensor
real(pReal), dimension(6) :: myTstar_v ! local copy of the 2nd Piola-Kirchhoff stress vector
real(pReal), dimension(constitutive_maxSizeState) :: myState ! local copy of the state
real(pReal), dimension(constitutive_maxSizeState) :: myState, & ! local copy of the state
myDotState ! local copy of dotState
integer(pInt) NiterationCrystallite, & ! number of iterations in crystallite loop
NiterationState ! number of iterations in state loop
integer(pInt) e, & ! element index
@ -267,8 +284,11 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
k, &
l, &
myNgrains, &
mySizeState
logical onTrack, & ! flag indicating wether we are still on track
mySizeState, &
mySizeDotState
logical onTrack, & ! flag indicating whether we are still on track
temperatureConverged, & ! flag indicating if temperature converged
stateConverged, & ! flag indicating if state converged
converged ! flag indicating if iteration converged
@ -315,8 +335,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
do while (any(crystallite_subStep(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) > subStepMin)) ! cutback loop for crystallites
NiterationCrystallite = NiterationCrystallite + 1
if (any(.not. crystallite_converged(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) .and. & ! any non-converged grain
.not. crystallite_localConstitution(:,:,FEsolving_execELem(1):FEsolving_execElem(2))) ) & ! has non-local constitution?
crystallite_converged(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) = &
@ -328,7 +346,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains
debugger = (e == 1 .and. i == 1 .and. g == 1)
if (crystallite_converged(g,i,e)) then
if (debugger) then
!$OMP CRITICAL (write2out)
@ -347,11 +364,17 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
crystallite_subLp0(:,:,g,i,e) = crystallite_Lp(:,:,g,i,e) ! ...plastic velocity gradient
constitutive_subState0(g,i,e)%p = constitutive_state(g,i,e)%p ! ...microstructure
crystallite_subTstar0_v(:,g,i,e) = crystallite_Tstar_v(:,g,i,e) ! ...2nd PK stress
else ! already at final time
!$OMP CRITICAL (distributionCrystallite)
debug_CrystalliteLoopDistribution(min(nCryst+1,NiterationCrystallite)) = &
debug_CrystalliteLoopDistribution(min(nCryst+1,NiterationCrystallite)) + 1
!$OMPEND CRITICAL (distributionCrystallite)
endif
else
crystallite_subStep(g,i,e) = 0.5_pReal * crystallite_subStep(g,i,e) ! cut step in half and restore...
crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) ! ...temperature
crystallite_Fp(:,:,g,i,e) = crystallite_subFp0(:,:,g,i,e) ! ...plastic def grad
crystallite_invFp(:,:,g,i,e) = math_inv3x3(crystallite_Fp(:,:,g,i,e))
crystallite_Lp(:,:,g,i,e) = crystallite_subLp0(:,:,g,i,e) ! ...plastic velocity grad
constitutive_state(g,i,e)%p = constitutive_subState0(g,i,e)%p ! ...microstructure
crystallite_Tstar_v(:,g,i,e) = crystallite_subTstar0_v(:,g,i,e) ! ...2nd PK stress
@ -383,36 +406,54 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
enddo
!$OMPEND PARALLEL DO
crystallite_nonfinished = ( crystallite_requested &
.and. crystallite_onTrack &
.and. .not. crystallite_converged)
! --+>> preguess for state <<+--
!
! incrementing by crystallite_subdt
! based on constitutive_subState0
! results in constitutive_state
! first loop for collection of state evolution based on old state
! second loop for updating to new state
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains
if ( crystallite_requested(g,i,e) &
.and. crystallite_onTrack(g,i,e) &
.and. .not. crystallite_converged(g,i,e)) then ! all undone crystallites
crystallite_converged(g,i,e) = crystallite_updateState(g,i,e) ! use former evolution rate
crystallite_converged(g,i,e) = crystallite_updateTemperature(g,i,e) ! use former evolution rate
constitutive_dotState(g,i,e)%p = 0.0_pReal ! zero out dotState
enddo; enddo; enddo
!$OMPEND PARALLEL DO
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains
if (crystallite_nonfinished(g,i,e)) & ! all undone crystallites
call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_Fp(:,:,g,i,e), &
crystallite_invFp(:,:,g,i,e), crystallite_Temperature(g,i,e), g, i, e)
enddo; enddo; enddo
!$OMPEND PARALLEL DO
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains
if (crystallite_nonfinished(g,i,e)) then ! all undone crystallites
crystallite_stateConverged(g,i,e) = crystallite_updateState(g,i,e) ! update state
crystallite_temperatureConverged(g,i,e) = crystallite_updateTemperature(g,i,e) ! update temperature
crystallite_converged(g,i,e) = .false. ! force at least one iteration step even if state already converged
endif
enddo
enddo
enddo
enddo; enddo; enddo
!$OMPEND PARALLEL DO
! --+>> state loop <<+--
NiterationState = 0_pInt
do while ( any( crystallite_requested(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) &
.and. crystallite_onTrack(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) &
.and. .not. crystallite_converged(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) ) &
do while ( any(crystallite_nonfinished(:,:,FEsolving_execELem(1):FEsolving_execElem(2))) &
.and. NiterationState < nState) ! convergence loop for crystallite
NiterationState = NiterationState + 1_pInt
@ -430,41 +471,57 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains
debugger = (e == 1 .and. i == 1 .and. g == 1)
if ( crystallite_requested(g,i,e) &
.and. crystallite_onTrack(g,i,e) &
.and. .not. crystallite_converged(g,i,e) ) & ! all undone crystallites
if (crystallite_nonfinished(g,i,e)) & ! all undone crystallites
crystallite_onTrack(g,i,e) = crystallite_integrateStress(g,i,e)
enddo
enddo
enddo
!$OMPEND PARALLEL DO
crystallite_nonfinished = crystallite_nonfinished .and. crystallite_onTrack
! --+>> state integration <<+--
!
! incrementing by crystallite_subdt
! based on constitutive_subState0
! results in constitutive_state
! first loop for collection of state evolution based on old state
! second loop for updating to new state
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains
constitutive_dotState(g,i,e)%p = 0.0_pReal ! zero out dotState
enddo; enddo; enddo
!$OMPEND PARALLEL DO
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains
if (crystallite_nonfinished(g,i,e)) & ! all undone crystallites
call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_Fp(:,:,g,i,e), &
crystallite_invFp(:,:,g,i,e), crystallite_Temperature(g,i,e), g, i, e)
enddo; enddo; enddo
!$OMPEND PARALLEL DO
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains
debugger = (e == 1 .and. i == 1 .and. g == 1)
if ( crystallite_requested(g,i,e) &
.and. crystallite_onTrack(g,i,e) &
.and. .not. crystallite_converged(g,i,e)) then ! all undone crystallites
crystallite_converged(g,i,e) = crystallite_updateState(g,i,e) .and. &
crystallite_updateTemperature(g,i,e)
if (crystallite_nonfinished(g,i,e)) then ! all undone crystallites
crystallite_stateConverged(g,i,e) = crystallite_updateState(g,i,e) ! update state
crystallite_temperatureConverged(g,i,e) = crystallite_updateTemperature(g,i,e) ! update temperature
crystallite_converged(g,i,e) = crystallite_stateConverged(g,i,e) .and. crystallite_temperatureConverged(g,i,e)
if (debugger) write (6,*) g,i,e,'converged after updState',crystallite_converged(g,i,e)
if (crystallite_converged(g,i,e)) then
!$OMP CRITICAL (distributionState)
debug_StateLoopDistribution(NiterationState) = debug_StateLoopDistribution(NiterationState) + 1
debug_CrystalliteStateLoopDistribution(NiterationState) = &
debug_CrystalliteStateLoopDistribution(NiterationState) + 1
!$OMPEND CRITICAL (distributionState)
!$OMP CRITICAL (distributionCrystallite)
debug_CrystalliteLoopDistribution(NiterationCrystallite) = &
debug_CrystalliteLoopDistribution(NiterationCrystallite) + 1
!$OMPEND CRITICAL (distributionCrystallite)
endif
endif
enddo
@ -472,8 +529,12 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
enddo
!$OMPEND PARALLEL DO
crystallite_nonfinished = crystallite_nonfinished .and. .not. crystallite_converged
enddo ! crystallite convergence loop
NiterationCrystallite = NiterationCrystallite + 1
enddo ! cutback loop
@ -511,10 +572,12 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
! debugger = (g == 1 .and. i == 1 .and. e == 1)
if (crystallite_converged(g,i,e)) then ! grain converged in above iteration
mySizeState = constitutive_sizeState(g,i,e) ! number of state variables for this grain
myState(1:mySizeState) = constitutive_state(g,i,e)%p ! remember unperturbed, converged state...
mySizeDotState = constitutive_sizeDotState(g,i,e) ! number of dotStates for this grain
myState(1:mySizeState) = constitutive_state(g,i,e)%p ! remember unperturbed, converged state, ...
myDotState(1:mySizeDotState) = constitutive_dotState(g,i,e)%p ! ... dotStates, ...
myTemperature = crystallite_Temperature(g,i,e) ! ... Temperature, ...
myF = crystallite_subF(:,:,g,i,e) ! ... and kinematics
myFp = crystallite_Fp(:,:,g,i,e)
myTemperature = crystallite_Temperature(g,i,e)
myFe = crystallite_Fe(:,:,g,i,e)
myLp = crystallite_Lp(:,:,g,i,e)
myTstar_v = crystallite_Tstar_v(:,g,i,e)
@ -548,8 +611,11 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
do while(.not. converged .and. onTrack .and. NiterationState < nState) ! keep cycling until done (potentially non-converged)
NiterationState = NiterationState + 1_pInt
onTrack = crystallite_integrateStress(g,i,e) ! stress of perturbed situation (overwrites _P,_Tstar_v,_Fp,_Lp,_Fe)
if (onTrack) converged = crystallite_updateState(g,i,e).AND.& ! update state
crystallite_updateTemperature(g,i,e) ! update temperature
if (onTrack) then
stateConverged = crystallite_updateState(g,i,e) ! update state
temperatureConverged = crystallite_updateTemperature(g,i,e) ! update temperature
converged = stateConverged .and. temperatureConverged
endif
if (debugger) then
!$OMP CRITICAL (write2out)
write (6,*) '-------------'
@ -563,8 +629,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
enddo
if (converged) & ! converged state warrants stiffness update
crystallite_dPdF(:,:,k,l,g,i,e) = (crystallite_P(:,:,g,i,e) - myP)/pert_Fg ! tangent dP_ij/dFg_kl
constitutive_state(g,i,e)%p = myState ! restore unperturbed, converged state...
crystallite_Temperature(g,i,e)= myTemperature ! ... temperature
constitutive_state(g,i,e)%p = myState ! restore unperturbed, converged state, ...
constitutive_dotState(g,i,e)%p = myDotState ! ... dotState, ...
crystallite_Temperature(g,i,e) = myTemperature ! ... temperature, ...
crystallite_Fp(:,:,g,i,e) = myFp ! ... and kinematics
crystallite_Fe(:,:,g,i,e) = myFe
crystallite_Lp(:,:,g,i,e) = myLp
@ -573,6 +640,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
!$OMP CRITICAL (out)
debug_StiffnessStateLoopDistribution(NiterationState) = &
debug_StiffnessstateLoopDistribution(NiterationState) + 1
if (nState < NiterationState) write(6,*) 'ohh shit!! stiffenss state loop debugging exceeded',NiterationState
!$OMPEND CRITICAL (out)
enddo
enddo
@ -633,9 +701,11 @@ endsubroutine
mySize = constitutive_sizeDotState(g,i,e)
! calculate the residuum
if (any(abs(constitutive_dotState(g,i,e)%p) > 1e10_pReal)) &
write(6,*) 'dotState: crap found at',g,i,e,constitutive_dotState(g,i,e)%p
call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
residuum = constitutive_state(g,i,e)%p(1:mySize) - constitutive_subState0(g,i,e)%p(1:mySize) - &
crystallite_subdt(g,i,e) * constitutive_dotState(crystallite_Tstar_v(:,g,i,e),crystallite_Temperature(g,i,e),g,i,e)
crystallite_subdt(g,i,e) * constitutive_dotState(g,i,e)%p(1:mySize)
call system_clock(count=tock,count_rate=tickrate,count_max=maxticks)
debug_cumDotStateCalls = debug_cumDotStateCalls + 1_pInt
debug_cumDotStateTicks = debug_cumDotStateTicks + tock-tick
@ -702,7 +772,8 @@ endsubroutine
! calculate the residuum
call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
residuum = crystallite_Temperature(g,i,e) - crystallite_subTemperature0(g,i,e) - &
crystallite_subdt(g,i,e) * constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e),crystallite_Temperature(g,i,e),g,i,e)
crystallite_subdt(g,i,e) * &
constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e),crystallite_Temperature(g,i,e),g,i,e)
call system_clock(count=tock,count_rate=tickrate,count_max=maxticks)
debug_cumDotTemperatureCalls = debug_cumDotTemperatureCalls + 1_pInt
debug_cumDotTemperatureTicks = debug_cumDotTemperatureTicks + tock-tick
@ -844,7 +915,7 @@ endsubroutine
A = math_mul33x33(transpose(invFp_current), math_mul33x33(transpose(Fg_new),math_mul33x33(Fg_new,invFp_current)))
! update microstructure
call constitutive_microstructure(crystallite_Temperature(g,i,e),g,i,e)
call constitutive_microstructure(crystallite_subTemperature0(g,i,e), crystallite_subFp0, g, i, e)
! get elasticity tensor
C_66 = constitutive_homogenizedC(g,i,e)
@ -877,7 +948,7 @@ LpLoop: do
! calculate plastic velocity gradient and its tangent according to constitutive law
call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
call constitutive_LpAndItsTangent(Lp_constitutive,dLp_constitutive,Tstar_v,crystallite_Temperature(g,i,e),g,i,e)
call constitutive_LpAndItsTangent(Lp_constitutive, dLp_constitutive, Tstar_v, crystallite_Temperature(g,i,e), g, i, e)
call system_clock(count=tock,count_rate=tickrate,count_max=maxticks)
debug_cumLpCalls = debug_cumLpCalls + 1_pInt
debug_cumLpTicks = debug_cumLpTicks + tock-tick
@ -934,6 +1005,10 @@ LpLoop: do
if (debugger) then
!$OMP CRITICAL (write2out)
write(6,*) '::: integrateStress failed on dR/dLp inversion at iteration', NiterationStress
write(6,'(a9,3(i3,x),/,9(9(e12.2,x)/))') 'dTdLp at ',g,i,e,dTdLp
write(6,'(a20,3(i3,x),/,9(9(e12.2,x)/))') 'dLp_constitutive at ',g,i,e,dLp_constitutive
write(6,'(a9,3(i3,x),/,9(9(f12.7,x)/))') 'dRdLp at ',g,i,e,dRdLp
write(6,'(a11,3(i3,x),/,3(3(f12.7,x)/))') 'Lpguess at ',g,i,e,Lpguess
!$OMPEND CRITICAL (write2out)
endif
return
@ -979,6 +1054,7 @@ LpLoop: do
crystallite_Tstar_v(:,g,i,e) = Tstar_v
crystallite_Fp(:,:,g,i,e) = Fp_new
crystallite_Fe(:,:,g,i,e) = Fe_new
crystallite_invFp(:,:,g,i,e) = invFp_new
! set return flag to true
crystallite_integrateStress = .true.
@ -993,6 +1069,7 @@ LpLoop: do
!$OMP CRITICAL (distributionStress)
debug_StressLoopDistribution(NiterationStress) = debug_StressLoopDistribution(NiterationStress) + 1
if (nStress < NiterationStress) write(6,*) 'ohh shit!! debug loop of stress exceeded',NiterationStress
!$OMPEND CRITICAL (distributionStress)
return

80
code/debug.f90
View File

@ -6,10 +6,11 @@
implicit none
integer(pInt), dimension(:), allocatable :: debug_StressLoopDistribution
integer(pInt), dimension(:), allocatable :: debug_StateLoopDistribution
integer(pInt), dimension(:), allocatable :: debug_CrystalliteStateLoopDistribution
integer(pInt), dimension(:), allocatable :: debug_StiffnessStateLoopDistribution
integer(pInt), dimension(:), allocatable :: debug_CrystalliteLoopDistribution
integer(pInt), dimension(:), allocatable :: debug_MaterialpointLoopDistribution ! added <<<updated 31.07.2009>>>
integer(pInt), dimension(:), allocatable :: debug_MaterialpointStateLoopDistribution
integer(pInt), dimension(:), allocatable :: debug_MaterialpointLoopDistribution
integer(pLongInt) :: debug_cumLpTicks = 0_pInt
integer(pLongInt) :: debug_cumDotStateTicks = 0_pInt
integer(pLongInt) :: debug_cumDotTemperatureTicks = 0_pInt
@ -27,7 +28,8 @@ subroutine debug_init()
use numerics, only: nStress, &
nState, &
nCryst, &
nHomog ! added <<<updated 31.07.2009>>>
nMPstate, &
nHomog
implicit none
write(6,*)
@ -35,10 +37,11 @@ subroutine debug_init()
write(6,*)
allocate(debug_StressLoopDistribution(nStress)) ; debug_StressLoopDistribution = 0_pInt
allocate(debug_StateLoopDistribution(nState)) ; debug_StateLoopDistribution = 0_pInt
allocate(debug_CrystalliteStateLoopDistribution(nState)); debug_CrystalliteStateLoopDistribution = 0_pInt
allocate(debug_StiffnessStateLoopDistribution(nState)) ; debug_StiffnessStateLoopDistribution = 0_pInt
allocate(debug_CrystalliteLoopDistribution(nCryst)) ; debug_CrystalliteLoopDistribution = 0_pInt
allocate(debug_MaterialpointLoopDistribution(nhomog)) ; debug_MaterialpointLoopDistribution = 0_pInt ! added <<<updated 31.07.2009>>>
allocate(debug_CrystalliteLoopDistribution(nCryst+1)) ; debug_CrystalliteLoopDistribution = 0_pInt
allocate(debug_MaterialpointStateLoopDistribution(nMPstate)) ; debug_MaterialpointStateLoopDistribution = 0_pInt
allocate(debug_MaterialpointLoopDistribution(nHomog+1)) ; debug_MaterialpointLoopDistribution = 0_pInt
endsubroutine
!********************************************************************
@ -50,10 +53,11 @@ subroutine debug_reset()
implicit none
debug_StressLoopDistribution = 0_pInt ! initialize debugging data
debug_StateLoopDistribution = 0_pInt
debug_CrystalliteStateLoopDistribution = 0_pInt
debug_StiffnessStateLoopDistribution = 0_pInt
debug_CrystalliteLoopDistribution = 0_pInt
debug_MaterialpointLoopDistribution = 0_pInt ! added <<<updated 31.07.2009>>>
debug_MaterialpointStateLoopDistribution = 0_pInt
debug_MaterialpointLoopDistribution = 0_pInt
debug_cumLpTicks = 0_pInt
debug_cumDotStateTicks = 0_pInt
debug_cumDotTemperatureTicks = 0_pInt
@ -72,7 +76,8 @@ endsubroutine
use numerics, only: nStress, &
nState, &
nCryst, &
nHomog ! added <<<updated 31.07.2009>>>
nMPstate, &
nHomog
implicit none
integer(pInt) i,integral
@ -111,21 +116,21 @@ endsubroutine
do i=1,nStress
if (debug_StressLoopDistribution(i) /= 0) then
integral = integral + i*debug_StressLoopDistribution(i)
write(6,'(i25,i10)') i,debug_StressLoopDistribution(i)
write(6,'(i25,x,i10)') i,debug_StressLoopDistribution(i)
endif
enddo
write(6,'(a15,i10,i10)') ' total',integral,sum(debug_StressLoopDistribution)
write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_StressLoopDistribution)
integral = 0_pInt
write(6,*)
write(6,*) 'distribution_StateLoop :'
write(6,*) 'distribution_CrystalliteStateLoop :'
do i=1,nState
if (debug_StateLoopDistribution(i) /= 0) then
integral = integral + i*debug_StateLoopDistribution(i)
write(6,'(i25,i10)') i,debug_StateLoopDistribution(i)
if (debug_CrystalliteStateLoopDistribution(i) /= 0) then
integral = integral + i*debug_CrystalliteStateLoopDistribution(i)
write(6,'(i25,x,i10)') i,debug_CrystalliteStateLoopDistribution(i)
endif
enddo
write(6,'(a15,i10,i10)') ' total',integral,sum(debug_StateLoopDistribution)
write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_CrystalliteStateLoopDistribution)
integral = 0_pInt
write(6,*)
@ -133,36 +138,57 @@ endsubroutine
do i=1,nState
if (debug_StiffnessStateLoopDistribution(i) /= 0) then
integral = integral + i*debug_StiffnessStateLoopDistribution(i)
write(6,'(i25,i10)') i,debug_StiffnessStateLoopDistribution(i)
write(6,'(i25,x,i10)') i,debug_StiffnessStateLoopDistribution(i)
endif
enddo
write(6,'(a15,i10,i10)') ' total',integral,sum(debug_StiffnessStateLoopDistribution)
write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_StiffnessStateLoopDistribution)
integral = 0_pInt
write(6,*)
write(6,*) 'distribution_CrystalliteLoop :'
do i=1,nCryst
do i=1,nCryst+1
if (debug_CrystalliteLoopDistribution(i) /= 0) then
integral = integral + i*debug_CrystalliteLoopDistribution(i)
write(6,'(i25,i10)') i,debug_CrystalliteLoopDistribution(i)
if (i <= nCryst) then
write(6,'(i25,x,i10)') i,debug_CrystalliteLoopDistribution(i)
else
write(6,'(i25,a1,i10)') i-1,'+',debug_CrystalliteLoopDistribution(i)
endif
endif
enddo
write(6,'(a15,i10,i10)') ' total',integral,sum(debug_CrystalliteLoopDistribution)
write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_CrystalliteLoopDistribution)
write(6,*)
!* Material point loop counter <<<updated 31.07.2009>>>
integral = 0_pInt
write(6,*)
write(6,*) 'distribution_MaterialpointLoop :'
do i=1,nCryst
if (debug_MaterialpointLoopDistribution(i) /= 0) then
integral = integral + i*debug_MaterialpointLoopDistribution(i)
write(6,'(i25,i10)') i,debug_MaterialpointLoopDistribution(i)
write(6,*) 'distribution_MaterialpointStateLoop :'
do i=1,nMPstate
if (debug_MaterialpointStateLoopDistribution(i) /= 0) then
integral = integral + i*debug_MaterialpointStateLoopDistribution(i)
write(6,'(i25,x,i10)') i,debug_MaterialpointStateLoopDistribution(i)
endif
enddo
write(6,'(a15,i10,i10)') ' total',integral,sum(debug_MaterialpointLoopDistribution)
write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_MaterialpointStateLoopDistribution)
write(6,*)
integral = 0_pInt
write(6,*)
write(6,*) 'distribution_MaterialpointLoop :'
do i=1,nHomog+1
if (debug_MaterialpointLoopDistribution(i) /= 0) then
integral = integral + i*debug_MaterialpointLoopDistribution(i)
if (i <= nHomog) then
write(6,'(i25,x,i10)') i,debug_MaterialpointLoopDistribution(i)
else
write(6,'(i25,a1,i10)') i-1,'+',debug_MaterialpointLoopDistribution(i)
endif
endif
enddo
write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_MaterialpointLoopDistribution)
write(6,*)
endsubroutine
END MODULE debug

55
code/homogenization.f90
View File

@ -185,9 +185,11 @@ subroutine materialpoint_stressAndItsTangent(&
use prec, only: pInt, &
pReal
use numerics, only: subStepMin, &
nHomog
nHomog, &
nMPstate
use FEsolving, only: FEsolving_execElem, &
FEsolving_execIP
FEsolving_execIP, &
terminallyIll
use mesh, only: mesh_element
use material, only: homogenization_Ngrains
use constitutive, only: constitutive_state0, &
@ -209,14 +211,16 @@ subroutine materialpoint_stressAndItsTangent(&
crystallite_partionedTstar0_v, &
crystallite_dt, &
crystallite_requested, &
crystallite_converged, &
crystallite_stressAndItsTangent
use debug, only: debug_MaterialpointLoopdistribution ! added <<<updated 31.07.2009>>>
use debug, only: debug_MaterialpointLoopDistribution, &
debug_MaterialpointStateLoopDistribution
implicit none
real(pReal), intent(in) :: dt
logical, intent(in) :: updateJaco
integer(pInt) homogenization_Niteration
integer(pInt) NiterationHomog,NiterationMPstate
integer(pInt) g,i,e,myNgrains
! ------ initialize to starting condition ------
@ -255,6 +259,7 @@ subroutine materialpoint_stressAndItsTangent(&
enddo
!$OMP END PARALLEL DO
NiterationHomog = 0_pInt
! ------ cutback loop ------
@ -285,7 +290,11 @@ subroutine materialpoint_stressAndItsTangent(&
if (homogenization_sizeState(i,e) > 0_pInt) &
homogenization_subState0(i,e)%p = homogenization_state(i,e)%p ! ...internal state of homog scheme
materialpoint_subF0(:,:,i,e) = materialpoint_subF(:,:,i,e) ! ...def grad
else ! already at final time
!$OMP CRITICAL (distributionHomog)
debug_MaterialpointLoopDistribution(min(nHomog+1,NiterationHomog)) = &
debug_MaterialpointLoopDistribution(min(nHomog+1,NiterationHomog)) + 1
!$OMPEND CRITICAL (distributionHomog)
endif
! materialpoint didn't converge, so we need a cutback here
@ -316,19 +325,19 @@ subroutine materialpoint_stressAndItsTangent(&
!$OMP END PARALLEL DO
!* Checks for cutback/substepping loops: added <<<updated 31.07.2009>>>
write (6,'(a,/,8(L,x))') 'MP exceeds substep min',materialpoint_subStep(:,FEsolving_execELem(1):FEsolving_execElem(2)) > subStepMin
write (6,'(a,/,8(L,x))') 'MP requested',materialpoint_requested(:,FEsolving_execELem(1):FEsolving_execElem(2))
write (6,'(a,/,8(f6.4,x))') 'MP subFrac',materialpoint_subFrac(:,FEsolving_execELem(1):FEsolving_execElem(2))
write (6,'(a,/,8(f6.4,x))') 'MP subStep',materialpoint_subStep(:,FEsolving_execELem(1):FEsolving_execElem(2))
! write (6,'(a,/,8(L,x))') 'MP exceeds substep min',materialpoint_subStep(:,FEsolving_execELem(1):FEsolving_execElem(2)) > subStepMin
! write (6,'(a,/,8(L,x))') 'MP requested',materialpoint_requested(:,FEsolving_execELem(1):FEsolving_execElem(2))
! write (6,'(a,/,8(f6.4,x))') 'MP subFrac',materialpoint_subFrac(:,FEsolving_execELem(1):FEsolving_execElem(2))
! write (6,'(a,/,8(f6.4,x))') 'MP subStep',materialpoint_subStep(:,FEsolving_execELem(1):FEsolving_execElem(2))
! ------ convergence loop material point homogenization ------
homogenization_Niteration = 0_pInt
NiterationMPstate = 0_pInt
do while (any( materialpoint_requested(:,FEsolving_execELem(1):FEsolving_execElem(2)) &
.and. .not. materialpoint_doneAndHappy(1,:,FEsolving_execELem(1):FEsolving_execElem(2)) &
) .and. homogenization_Niteration < nHomog) ! convergence loop for materialpoint
homogenization_Niteration = homogenization_Niteration + 1
) .and. NiterationMPstate < nMPstate) ! convergence loop for materialpoint
NiterationMPstate = NiterationMPstate + 1
! --+>> deformation partitioning <<+--
!
@ -367,11 +376,15 @@ subroutine materialpoint_stressAndItsTangent(&
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
if ( materialpoint_requested(i,e) .and. &
.not. materialpoint_doneAndHappy(1,i,e)) then
if (.not. all(crystallite_converged(:,i,e))) then
materialpoint_doneAndHappy(:,i,e) = (/.true.,.false./)
else
materialpoint_doneAndHappy(:,i,e) = homogenization_updateState(i,e)
endif
materialpoint_converged(i,e) = all(materialpoint_doneAndHappy(:,i,e)) ! converged if done and happy
if (materialpoint_converged(i,e)) & ! added <<<updated 31.07.2009>>>
debug_MaterialpointLoopdistribution(homogenization_Niteration) = &
debug_MaterialpointLoopdistribution(homogenization_Niteration) + 1
debug_MaterialpointStateLoopdistribution(NiterationMPstate) = &
debug_MaterialpointStateLoopdistribution(NiterationMPstate) + 1
endif
enddo
enddo
@ -379,18 +392,26 @@ subroutine materialpoint_stressAndItsTangent(&
enddo ! homogenization convergence loop
NiterationHomog = NiterationHomog +1_pInt
enddo ! cutback loop
! check for non-performer: any(.not. converged)
! replace with elastic response ?
! replace everybody with odd response ?
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
elementLoop: do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
if (materialpoint_converged(i,e)) then
call homogenization_averageStressAndItsTangent(i,e)
call homogenization_averageTemperature(i,e)
else
terminallyIll = .true.
exit elementLoop
endif
enddo
enddo
enddo elementLoop
!$OMP END PARALLEL DO
write (6,*) 'Material Point finished'

3
code/homogenization_RGC.f90
View File

@ -64,8 +64,7 @@ subroutine homogenization_RGC_init(&
allocate(homogenization_RGC_Ngrains(3,maxNinstance)); homogenization_RGC_Ngrains = 0_pInt
allocate(homogenization_RGC_ciAlpha(3,maxNinstance)); homogenization_RGC_ciAlpha = 0.0_pReal
allocate(homogenization_RGC_xiAlpha(3,maxNinstance)); homogenization_RGC_xiAlpha = 0.0_pReal
allocate(homogenization_RGC_output(maxval(homogenization_Noutput),maxNinstance)); &
homogenization_RGC_output = ''
allocate(homogenization_RGC_output(maxval(homogenization_Noutput),maxNinstance)); homogenization_RGC_output = ''
rewind(file)
line = ''

97
code/material.config
View File

@ -6,56 +6,23 @@
type isostrain
Ngrains 1
[RGC]
type RGC
Ngrains 8
[Taylor2]
type isostrain
Ngrains 2
[Taylor4]
type isostrain
Ngrains 4
[Taylor10]
type isostrain
Ngrains 10
[Taylor100]
type isostrain
Ngrains 100
#####################
<microstructure>
#####################
[TWIPSteel_Poly]
(constituent) phase 5 texture 1 fraction 1.0
[Aluminum_Poly]
(constituent) phase 3 texture 1 fraction 1.0
[TWIPsteel_001]
(constituent) phase 5 texture 2 fraction 1.0
[Aluminum_001]
(constituent) phase 3 texture 2 fraction 1.0
[TWIPsteel_011]
(constituent) phase 5 texture 3 fraction 1.0
[TWIPsteel_111]
(constituent) phase 5 texture 4 fraction 1.0
[TWIPsteel_123]
(constituent) phase 5 texture 5 fraction 1.0
[Alu_Polycrystal]
[Aluminum_j2]
(constituent) phase 1 texture 1 fraction 1.0
[Alu_001]
(constituent) phase 1 texture 2 fraction 1.0
[Alu_011]
(constituent) phase 1 texture 3 fraction 1.0
[Alu_111]
(constituent) phase 1 texture 4 fraction 1.0
#####################
<phase>
@ -101,11 +68,11 @@ c44 28.34e9
gdot0_slip 0.001
n_slip 20
s0_slip 31e6 0 0 0 # per family
ssat_slip 63e6 0 0 0 # per family
tau0_slip 31e6 # per family
tausat_slip 63e6 # per family
gdot0_twin 0.001
n_twin 20
s0_twin 31e6 0 0 0 # per family
tau0_twin 31e6 # per family
s_pr 0 # push-up factor for slip saturation due to twinning
twin_b 0
twin_c 0
@ -120,43 +87,27 @@ interaction_sliptwin 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
interaction_twinslip 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
interaction_twintwin 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[Aluminum_nonlocal]
[TWIP steel FeMnC]
constitution nonlocal
/nonlocal/
constitution dislobased
(output) dislocation_density
(output) state_slip
(output) rateofshear_slip
(output) mfp_slip
(output) thresholdstress_slip
(output) state_twin
(output) rateofshear_twin
(output) mfp_twin
(output) thresholdstress_twin
C11 175.0e9 # elastic constants in Pa
C12 115.0e9
C44 135.0e9
lattice_structure fcc
Nslip 12
Ntwin 12
Nslip 12 0 0 0 # per family
### dislocation density-based constitutive parameters ###
burgers 2.56e-10 # Burgers vector [m]
Qedge 5.5e-19 # Activation energy for dislocation glide [J/K] (0.5*G*b^3)
grainsize 2.0e-5 # Average grain size [m]
stacksize 5.0e-8 # Twin stack mean thickness [m]
fmax 1.0 # Maximum admissible twin volume fraction
Ndot0 0.0 # Number of potential twin source per volume per time [1/m³.s]
interaction_coefficients 1.0 2.2 3.0 1.6 3.8 4.5 # Dislocation interaction coefficients
rho0 6.0e12 # Initial dislocation density [m/m³]
Cmfpslip 2.0 # Adjustable parameter controlling dislocation mean free path
Cmfptwin 1.0 # Adjustable parameter controlling twin mean free path
Cthresholdslip 0.1 # Adjustable parameter controlling threshold stress for dislocation motion
Cthresholdtwin 1.0 # Adjustable parameter controlling threshold stress for deformation twinning
Cactivolume 1.0 # Adjustable parameter controlling activation volume
Cstorage 0.02 # Adjustable parameter controlling dislocation storage
Carecovery 0.0 # Adjustable parameter controlling athermal recovery
c11 106.75e9
c12 60.41e9
c44 28.34e9
burgers 2.56e-10 0 0 0 # Burgers vector in m
rhoEdgePos0 2.5e12 0 0 0 # Initial positive edge dislocation density in m/m**3
rhoEdgeNeg0 2.5e12 0 0 0 # Initial negative edge dislocation density in m/m**3
rhoScrewPos0 2.5e12 0 0 0 # Initial positive screw dislocation density in m/m**3
rhoScrewNeg0 2.5e12 0 0 0 # Initial negative screw dislocation density in m/m**3
v0 100 0 0 0 # initial dislocation velocity
interaction_SlipSlip 1.0 2.2 3.0 1.6 3.8 4.5 # Dislocation interaction coefficients
#####################
<texture>

100
code/math.f90
View File

@ -556,7 +556,6 @@ endsubroutine
!**************************************************************************
! matrix multiplication 3x3 = 1
!**************************************************************************
PURE FUNCTION math_mul3x3(A,B)
use prec, only: pReal, pInt
@ -578,7 +577,6 @@ endsubroutine
!**************************************************************************
! matrix multiplication 6x6 = 1
!**************************************************************************
PURE FUNCTION math_mul6x6(A,B)
use prec, only: pReal, pInt
@ -596,10 +594,10 @@ endsubroutine
ENDFUNCTION
!**************************************************************************
! matrix multiplication 33x33 = 3x3
!**************************************************************************
!**************************************************************************´
PURE FUNCTION math_mul33x33(A,B)
use prec, only: pReal, pInt
@ -636,27 +634,6 @@ endsubroutine
ENDFUNCTION
!**************************************************************************
! matrix multiplication 66x6 = 6
!**************************************************************************
PURE FUNCTION math_mul66x6(A,B)
use prec, only: pReal, pInt
implicit none
integer(pInt) i
real(pReal), dimension(6,6), intent(in) :: A
real(pReal), dimension(6), intent(in) :: B
real(pReal), dimension(6) :: math_mul66x6
forall (i=1:6) math_mul66x6(i) = &
A(i,1)*B(1) + A(i,2)*B(2) + A(i,3)*B(3) + &
A(i,4)*B(4) + A(i,5)*B(5) + A(i,6)*B(6)
return
ENDFUNCTION
!**************************************************************************
! matrix multiplication 99x99 = 9x9
!**************************************************************************
@ -680,6 +657,62 @@ endsubroutine
ENDFUNCTION
!**************************************************************************
! matrix multiplication 33x3 = 3
!**************************************************************************
PURE FUNCTION math_mul33x3(A,B)
use prec, only: pReal, pInt
implicit none
integer(pInt) i
real(pReal), dimension(3,3), intent(in) :: A
real(pReal), dimension(3), intent(in) :: B
real(pReal), dimension(3) :: math_mul33x3
forall (i=1:3) math_mul33x3(i) = A(i,1)*B(1) + A(i,2)*B(2) + A(i,3)*B(3)
return
ENDFUNCTION
!**************************************************************************
! matrix multiplication 66x6 = 6
!**************************************************************************
PURE FUNCTION math_mul66x6(A,B)
use prec, only: pReal, pInt
implicit none
integer(pInt) i
real(pReal), dimension(6,6), intent(in) :: A
real(pReal), dimension(6), intent(in) :: B
real(pReal), dimension(6) :: math_mul66x6
forall (i=1:6) math_mul66x6(i) = &
A(i,1)*B(1) + A(i,2)*B(2) + A(i,3)*B(3) + &
A(i,4)*B(4) + A(i,5)*B(5) + A(i,6)*B(6)
return
ENDFUNCTION
!**************************************************************************
! transposition of a 3x3 matrix
!**************************************************************************
function math_transpose3x3(A)
use prec, only: pReal,pInt
implicit none
real(pReal),dimension(3,3),intent(in) :: A
real(pReal),dimension(3,3) :: math_transpose3x3
math_transpose3x3 = reshape( (/A(1,1), A(1,2), A(1,3), A(2,1), A(2,2), A(2,3), A(3,1), A(3,2), A(3,3) /),(/3,3/) )
return
endfunction
!**************************************************************************
! Cramer inversion of 3x3 matrix (function)
@ -1025,6 +1058,23 @@ endsubroutine
ENDFUNCTION
!********************************************************************
! euclidic norm of a 3x1 vector
!********************************************************************
pure function math_norm3(v3)
use prec, only: pReal,pInt
implicit none
real(pReal), dimension(3), intent(in) :: v3
real(pReal) math_norm3
math_norm3 = sqrt(v3(1)*v3(1)+v3(2)*v3(2)+v3(3)*v3(3))
return
endfunction
!********************************************************************
! convert 3x3 matrix into vector 9x1
!********************************************************************

8
code/mesh.f90
View File

@ -46,13 +46,14 @@
integer(pInt), dimension(2) :: mesh_maxValStateVar = 0_pInt
character(len=64), dimension(:), allocatable :: mesh_nameElemSet
integer(pInt), dimension(:,:), allocatable :: mesh_mapElemSet
integer(pInt), dimension(:,:), allocatable, target :: mesh_mapFEtoCPelem,mesh_mapFEtoCPnode
integer(pInt), dimension(:,:), allocatable, target :: mesh_mapFEtoCPelem, mesh_mapFEtoCPnode
integer(pInt), dimension(:,:), allocatable :: mesh_element, mesh_sharedElem
integer(pInt), dimension(:,:,:,:), allocatable :: mesh_ipNeighborhood
real(pReal), dimension(:,:,:), allocatable :: mesh_subNodeCoord ! coordinates of subnodes per element
real(pReal), dimension(:,:), allocatable :: mesh_ipVolume ! volume associated with IP
real(pReal), dimension(:,:,:), allocatable :: mesh_ipArea ! area of interface to neighboring IP
real(pReal), dimension(:,:,:), allocatable :: mesh_ipArea, & ! area of interface to neighboring IP
mesh_ipCenterOfGravity ! center of gravity of IP
real(pReal), dimension(:,:,:,:), allocatable :: mesh_ipAreaNormal ! area normal of interface to neighboring IP
real(pReal), allocatable :: mesh_node (:,:)
@ -1578,6 +1579,8 @@ subroutine mesh_get_nodeElemDimensions (unit)
real(pReal), dimension(3) :: centerOfGravity
allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems)) ; mesh_ipVolume = 0.0_pReal
allocate(mesh_ipCenterOfGravity(3,mesh_maxNips,mesh_NcpElems)) ; mesh_ipCenterOfGravity = 0.0_pReal
do e = 1,mesh_NcpElems ! loop over cpElems
t = mesh_element(2,e) ! get elemType
do i = 1,FE_Nips(t) ! loop over IPs of elem
@ -1613,6 +1616,7 @@ subroutine mesh_get_nodeElemDimensions (unit)
mesh_ipVolume(i,e) = mesh_ipVolume(i,e) + sum(volume) ! add contribution from this interface
enddo
mesh_ipVolume(i,e) = mesh_ipVolume(i,e) / FE_NipFaceNodes ! renormalize with interfaceNodeNum due to loop over them
mesh_ipCenterOfGravity(:,i,e) = centerOfGravity
enddo
enddo
return

6
code/mpie_cpfem_marc.f90
View File

@ -46,6 +46,7 @@
include "constitutive_phenopowerlaw.f90" ! uses prec, math, IO, lattice, material, debug
include "constitutive_j2.f90" ! uses prec, math, IO, lattice, material, debug
include "constitutive_dislobased.f90" ! uses prec, math, IO, lattice, material, debug
include "constitutive_nonlocal.f90" ! uses prec, math, IO, lattice, material, debug
include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug
include "crystallite.f90" ! uses prec, math, IO, numerics
include "homogenization_isostrain.f90" ! uses prec, math, IO,
@ -166,7 +167,10 @@ subroutine hypela2(&
if (inc == 0) then
cycleCounter = 4
else
if (theCycle > ncycle .or. theInc /= inc) cycleCounter = 0 ! reset counter for each cutback or new inc
if (theCycle > ncycle .or. theInc /= inc) then
cycleCounter = 0 ! reset counter for each cutback or new inc
terminallyIll = .false.
endif
if (theCycle /= ncycle .or. theLovl /= lovl) then
cycleCounter = cycleCounter+1 ! ping pong
outdatedFFN1 = .false.

3
code/numerics.config
View File

@ -5,7 +5,8 @@ relevantStrain 1.0e-7 # strain increment considered signific
iJacoStiffness 1 # frequency of stiffness update
iJacoLpresiduum 1 # frequency of Jacobian update of residuum in Lp
pert_Fg 1.0e-6 # strain perturbation for FEM Jacobi
nHomog 10 # homogenization loop limit
nHomog 20 # homogenization loop limit
nMPstate 10 # material point state loop limit
nCryst 20 # crystallite loop limit (only for debugging info, real loop limit is "subStepMin")
nState 10 # state loop limit
nStress 40 # stress loop limit

12
code/numerics.f90
View File

@ -9,6 +9,7 @@ character(len=64), parameter :: numerics_configFile = 'numerics.config' ! name o
integer(pInt) iJacoStiffness, & ! frequency of stiffness update
iJacoLpresiduum, & ! frequency of Jacobian update of residuum in Lp
nHomog, & ! homogenization loop limit
nMPstate, & ! materialpoint state loop limit
nCryst, & ! crystallite loop limit (only for debugging info, real loop limit is "subStepMin")
nState, & ! state loop limit
nStress ! stress loop limit
@ -69,7 +70,8 @@ subroutine numerics_init()
iJacoStiffness = 1_pInt
iJacoLpresiduum = 1_pInt
pert_Fg = 1.0e-6_pReal
nHomog = 10_pInt
nHomog = 20_pInt
nMPstate = 10_pInt
nCryst = 20_pInt
nState = 10_pInt
nStress = 40_pInt
@ -111,6 +113,8 @@ subroutine numerics_init()
pert_Fg = IO_floatValue(line,positions,2)
case ('nhomog')
nHomog = IO_intValue(line,positions,2)
case ('nmpstate')
nMPstate = IO_intValue(line,positions,2)
case ('ncryst')
nCryst = IO_intValue(line,positions,2)
case ('nstate')
@ -160,6 +164,7 @@ subroutine numerics_init()
write(6,'(a24,x,i8)') 'iJacoLpresiduum: ',iJacoLpresiduum
write(6,'(a24,x,e8.1)') 'pert_Fg: ',pert_Fg
write(6,'(a24,x,i8)') 'nHomog: ',nHomog
write(6,'(a24,x,i8)') 'nMPstate: ',nMPstate
write(6,'(a24,x,i8)') 'nCryst: ',nCryst
write(6,'(a24,x,i8)') 'nState: ',nState
write(6,'(a24,x,i8)') 'nStress: ',nStress
@ -184,12 +189,13 @@ subroutine numerics_init()
if (iJacoLpresiduum < 1_pInt) call IO_error(262)
if (pert_Fg <= 0.0_pReal) call IO_error(263)
if (nHomog < 1_pInt) call IO_error(264)
if (nMPstate < 1_pInt) call IO_error(279) !! missing in IO !!
if (nCryst < 1_pInt) call IO_error(265)
if (nState < 1_pInt) call IO_error(266)
if (nStress < 1_pInt) call IO_error(267)
if (subStepMin <= 0.0_pReal) call IO_error(268)
if (rTol_crystalliteState <= 0.0_pReal) call IO_error(269)
if (rTol_crystalliteTemperature <= 0.0_pReal) call IO_error(276)
if (rTol_crystalliteTemperature <= 0.0_pReal) call IO_error(276) !! oops !!
if (rTol_crystalliteStress <= 0.0_pReal) call IO_error(270)
if (aTol_crystalliteStress <= 0.0_pReal) call IO_error(271)
@ -198,7 +204,7 @@ subroutine numerics_init()
if (relTol_RGC <= 0.0_pReal) call IO_error(273)
if (absMax_RGC <= 0.0_pReal) call IO_error(274)
if (relMax_RGC <= 0.0_pReal) call IO_error(275)
if (pPert_RGC <= 0.0_pReal) call IO_error(276)
if (pPert_RGC <= 0.0_pReal) call IO_error(276) !! oops !!
if (xSmoo_RGC <= 0.0_pReal) call IO_error(277)
endsubroutine