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

major restructuring of code. 

homogenization as well as constitutive are now free to choose. the runtime got somewhat longer (25% on simple tests) compared to a hardcoded isostrain homogenization. this might be a point of further optimization at a later stage...

please use homogenization_isostrain.f90 as starting point / example for future developments of homog-schemes.

the homogenization scheme now can additionally output certain results. hence, the userdata structure at each integration point now looks like this:
- sizeHomogPostResults
- block of that size containing homogPostResults
then for each grain:
- sizeGrainPostResults
- block of that size containing crystallitePostResults, which consist of:
+ phaseID
+ volFrac
+ Eulers (3)
+ any constitutive post results requested
This commit is contained in:
Philip Eisenlohr 2009-05-07 16:27:36 +00:00
parent 449b791463
commit 4f705f524f
14 changed files with 1913 additions and 906 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

774
trunk/CPFEM.f90
View File

@ -9,33 +9,10 @@
! **************************************************************** ! ****************************************************************
! *** General variables for the material behaviour calculation *** ! *** General variables for the material behaviour calculation ***
! **************************************************************** ! ****************************************************************
real(pReal), dimension (:,:), allocatable :: CPFEM_Temperature real(pReal), dimension (:,:,:), allocatable :: CPFEM_cs ! Cauchy stress
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_ffn_bar !average FFN per IP real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_dcsdE ! Cauchy stress tangent
real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_ffn !individual FFN per grain real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_dcsdE_knownGood ! known good tangent
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_ffn1_bar !average FFN1 per IP
real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_ffn1 !individual FFN1 per grain
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_PK1_bar !average PK1 per IP
real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_PK1 !individual PK1 per grain
real(pReal), dimension (:,:,:,:,:,:), allocatable :: CPFEM_dPdF_bar !average dPdF per IP
real(pReal), dimension (:,:,:,:,:,:), allocatable :: CPFEM_dPdF_bar_old !old average dPdF per IP
real(pReal), dimension (:,:,:,:,:,:,:),allocatable :: CPFEM_dPdF !individual dPdF per grain
real(pReal), dimension (:,:,:), allocatable :: CPFEM_stress_bar
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_jaco_bar
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_jaco_knownGood
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_results
real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_Lp_old
real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_Lp_new
real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_Fp_old
real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_Fp_new
real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_Fe_new
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_Tstar_v
logical, dimension (:,:,:), allocatable :: crystallite_converged !individual convergence flag per grain
integer(pInt), dimension(:,:), allocatable :: CPFEM_execution_IP
integer(pInt), dimension(2) :: CPFEM_execution_elem
integer(pInt) :: CPFEM_Nresults = 5_pInt ! phase, volfrac, three Euler angles
logical :: CPFEM_init_done = .false. ! remember whether init has been done already logical :: CPFEM_init_done = .false. ! remember whether init has been done already
logical :: CPFEM_calc_done = .false. ! remember whether first IP has already calced the results logical :: CPFEM_calc_done = .false. ! remember whether first IP has already calced the results
@ -47,89 +24,35 @@
!*** allocate the arrays defined in module CPFEM *** !*** allocate the arrays defined in module CPFEM ***
!*** and initialize them *** !*** and initialize them ***
!********************************************************* !*********************************************************
SUBROUTINE CPFEM_init(Temperature) SUBROUTINE CPFEM_init()
!
use prec use prec, only: pInt,pReal
use math, only: math_EulertoR, math_I3, math_identity2nd use FEsolving, only: parallelExecution,symmetricSolver,FEsolving_execElem,FEsolving_execIP
use FEsolving, only: parallelExecution use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips,FE_Nips
use mesh use material, only: homogenization_maxNgrains
use material use constitutive, only: constitutive_maxSizePostResults
use constitutive use crystallite, only: crystallite_Nresults
! use homogenization, only: homogenization_maxSizePostResults
implicit none implicit none
!
real(pReal) Temperature
integer(pInt) e,i,g integer(pInt) e,i,g
!
! *** mpie.marc parameters ***
allocate(CPFEM_Temperature(mesh_maxNips,mesh_NcpElems)) ; CPFEM_Temperature = Temperature
allocate(CPFEM_ffn_bar(3,3,mesh_maxNips,mesh_NcpElems))
forall(e=1:mesh_NcpElems,i=1:mesh_maxNips) CPFEM_ffn_bar(:,:,i,e) = math_I3
allocate(CPFEM_ffn(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
forall(g=1:homogenization_maxNgrains,e=1:mesh_NcpElems,i=1:mesh_maxNips) CPFEM_ffn(:,:,g,i,e) = math_I3
allocate(CPFEM_ffn1_bar(3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_ffn1_bar = CPFEM_ffn_bar
allocate(CPFEM_ffn1(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_ffn1 = CPFEM_ffn
allocate(CPFEM_PK1_bar(3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_PK1_bar = 0.0_pReal
allocate(CPFEM_PK1(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_PK1 = 0.0_pReal
allocate(CPFEM_dPdF_bar(3,3,3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dPdF_bar = 0.0_pReal
allocate(CPFEM_dPdF_bar_old(3,3,3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dPdF_bar_old = 0.0_pReal
allocate(CPFEM_dPdF(3,3,3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dPdF = 0.0_pReal
allocate(CPFEM_stress_bar(6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_stress_bar = 0.0_pReal
allocate(CPFEM_jaco_bar(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_jaco_bar = 0.0_pReal
allocate(CPFEM_jaco_knownGood(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_jaco_knownGood = 0.0_pReal
!
! *** User defined results ***
allocate(CPFEM_results(CPFEM_Nresults+constitutive_maxSizePostResults,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
CPFEM_results = 0.0_pReal
!
! *** Plastic velocity gradient ***
allocate(CPFEM_Lp_old(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Lp_old = 0.0_pReal
allocate(CPFEM_Lp_new(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Lp_new = 0.0_pReal
! *** Plastic deformation gradient at (t=t0) and (t=t1) *** allocate(CPFEM_cs(6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_cs = 0.0_pReal
allocate(CPFEM_Fp_new(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Fp_new = 0.0_pReal allocate(CPFEM_dcsdE(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dcsde = 0.0_pReal
allocate(CPFEM_Fp_old(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) allocate(CPFEM_dcsdE_knownGood(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dcsde_knownGood = 0.0_pReal
forall (e=1:mesh_NcpElems,i=1:mesh_maxNips,g=1:homogenization_maxNgrains) &
CPFEM_Fp_old(:,:,g,i,e) = math_EulerToR(material_EulerAngles(:,g,i,e)) ! plastic def gradient reflects init orientation
! *** Elastic deformation gradient at (t=t1) ***
allocate(CPFEM_Fe_new(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Fe_new = 0.0_pReal
! *** Stress vector at (t=t1) ***
allocate(CPFEM_Tstar_v(6,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Tstar_v = 0.0_pReal
!
allocate(crystallite_converged(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)); crystallite_converged = .false.
allocate(CPFEM_execution_IP(2,mesh_NcpElems)); CPFEM_execution_IP = 1_pInt
forall (e = 1:mesh_NcpElems) CPFEM_execution_IP(2,e) = FE_Nips(mesh_element(2,e))
CPFEM_execution_elem = (/1,mesh_NcpElems/)
! *** Output to MARC output file *** ! *** Output to MARC output file ***
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*) write(6,*)
write(6,*) 'CPFEM Initialization' write(6,*) '<<<+- cpfem init -+>>>'
write(6,*) write(6,*)
write(6,*) 'CPFEM_Temperature: ', shape(CPFEM_Temperature) write(6,'(a32,x,6(i5,x))') 'CPFEM_cs: ', shape(CPFEM_cs)
write(6,*) 'CPFEM_ffn_bar: ', shape(CPFEM_ffn_bar) write(6,'(a32,x,6(i5,x))') 'CPFEM_dcsde: ', shape(CPFEM_dcsde)
write(6,*) 'CPFEM_ffn: ', shape(CPFEM_ffn) write(6,'(a32,x,6(i5,x))') 'CPFEM_dcsde_knownGood: ', shape(CPFEM_dcsde_knownGood)
write(6,*) 'CPFEM_ffn1_bar: ', shape(CPFEM_ffn1_bar)
write(6,*) 'CPFEM_ffn1: ', shape(CPFEM_ffn1)
write(6,*) 'CPFEM_PK1_bar: ', shape(CPFEM_PK1_bar)
write(6,*) 'CPFEM_PK1: ', shape(CPFEM_PK1)
write(6,*) 'CPFEM_dPdF_bar: ', shape(CPFEM_dPdF_bar)
write(6,*) 'CPFEM_dPdF_bar_old: ', shape(CPFEM_dPdF_bar_old)
write(6,*) 'CPFEM_dPdF: ', shape(CPFEM_dPdF)
write(6,*) 'CPFEM_stress_bar: ', shape(CPFEM_stress_bar)
write(6,*) 'CPFEM_jaco_bar: ', shape(CPFEM_jaco_bar)
write(6,*) 'CPFEM_jaco_knownGood: ', shape(CPFEM_jaco_knownGood)
write(6,*) 'CPFEM_results: ', shape(CPFEM_results)
write(6,*) 'CPFEM_Lp_old: ', shape(CPFEM_Lp_old)
write(6,*) 'CPFEM_Lp_new: ', shape(CPFEM_Lp_new)
write(6,*) 'CPFEM_Fp_old: ', shape(CPFEM_Fp_old)
write(6,*) 'CPFEM_Fp_new: ', shape(CPFEM_Fp_new)
write(6,*) 'CPFEM_Fe_new: ', shape(CPFEM_Fe_new)
write(6,*) 'CPFEM_Tstar_v: ', shape(CPFEM_Tstar_v)
write(6,*) 'crystallite_converged:', shape(crystallite_converged)
write(6,*) write(6,*)
write(6,*) 'parallelExecution: ', parallelExecution write(6,*) 'parallelExecution: ', parallelExecution
write(6,*) 'symmetricSolver: ', symmetricSolver
call flush(6) call flush(6)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
return return
@ -153,7 +76,7 @@
! CPFEM_jaco jacobian in Mandel notation ! CPFEM_jaco jacobian in Mandel notation
! CPFEM_ngens size of stress strain law ! CPFEM_ngens size of stress strain law
!*********************************************************************** !***********************************************************************
SUBROUTINE CPFEM_general(CPFEM_mode, ffn, ffn1, Temperature, CPFEM_dt,& subroutine CPFEM_general(CPFEM_mode, ffn, ffn1, Temperature, CPFEM_dt,&
CPFEM_en, CPFEM_in, CPFEM_stress, CPFEM_updateJaco, CPFEM_jaco, CPFEM_ngens) CPFEM_en, CPFEM_in, CPFEM_stress, CPFEM_updateJaco, CPFEM_jaco, CPFEM_ngens)
! note: CPFEM_stress = Cauchy stress cs(6) and CPFEM_jaco = Consistent tangent dcs/de ! note: CPFEM_stress = Cauchy stress cs(6) and CPFEM_jaco = Consistent tangent dcs/de
! !
@ -161,15 +84,19 @@
use FEsolving use FEsolving
use debug use debug
use math use math
use mesh, only: mesh_init,mesh_FEasCP, mesh_NcpElems, mesh_maxNips, mesh_element use mesh, only: mesh_init,&
mesh_FEasCP,mesh_element,mesh_NcpElems,mesh_maxNips,FE_Nips
use lattice, only: lattice_init use lattice, only: lattice_init
use material use material, only: material_init, homogenization_maxNgrains
use constitutive, only: constitutive_init,constitutive_state_old,constitutive_state_new use constitutive, only: constitutive_init,&
constitutive_state0,constitutive_state
use crystallite
use homogenization
implicit none implicit none
!
integer(pInt) CPFEM_en, CPFEM_in, cp_en, CPFEM_ngens, i,j,k,l,m,n integer(pInt) CPFEM_en, CPFEM_in, cp_en, CPFEM_ngens, i,j,k,l,m,n
real(pReal), dimension (3,3) :: ffn,ffn1,Kirchhoff_bar real(pReal), dimension (3,3) :: ffn,ffn1,Kirchhoff
real(pReal), dimension (3,3,3,3) :: H_bar, H_bar_sym real(pReal), dimension (3,3,3,3) :: H, H_sym
real(pReal), dimension(CPFEM_ngens) :: CPFEM_stress real(pReal), dimension(CPFEM_ngens) :: CPFEM_stress
real(pReal), dimension(CPFEM_ngens,CPFEM_ngens) :: CPFEM_jaco real(pReal), dimension(CPFEM_ngens,CPFEM_ngens) :: CPFEM_jaco
real(pReal) Temperature,CPFEM_dt,J_inverse real(pReal) Temperature,CPFEM_dt,J_inverse
@ -179,627 +106,116 @@
! 4: recycling of former results (MARC speciality)& ! 4: recycling of former results (MARC speciality)&
! 5: record tangent from former converged inc& ! 5: record tangent from former converged inc&
! 6: restore tangent from former converged inc ! 6: restore tangent from former converged inc
integer(pInt) e
logical CPFEM_updateJaco logical CPFEM_updateJaco
!
if (.not. CPFEM_init_done) then ! initialization step (three dimensional stress state check missing?) if (.not. CPFEM_init_done) then ! initialization step (three dimensional stress state check missing?)
call math_init() call math_init()
call FE_init() call FE_init()
call mesh_init() call mesh_init()
FEsolving_execElem = (/1,mesh_NcpElems/)
allocate(FEsolving_execIP(2,mesh_NcpElems)); FEsolving_execIP = 1_pInt
forall (e = 1:mesh_NcpElems) FEsolving_execIP(2,e) = FE_Nips(mesh_element(2,e))
call lattice_init() call lattice_init()
call material_init() call material_init()
call constitutive_init() call constitutive_init()
write (6,*) 'call CPFEM init' call crystallite_init()
call CPFEM_init(Temperature) call homogenization_init()
call CPFEM_init()
CPFEM_init_done = .true. CPFEM_init_done = .true.
endif endif
!
if ((.not. parallelExecution) .and. (CPFEM_mode == 3)) CPFEM_mode = 2
!
cp_en = mesh_FEasCP('elem',CPFEM_en) cp_en = mesh_FEasCP('elem',CPFEM_en)
if (cp_en == 1 .and. CPFEM_in == 1) then if (cp_en == 1 .and. CPFEM_in == 1) then
write(6,*) '#####################################'
write(6,'(a10,1x,f8.4,1x,a10,1x,i4,1x,a10,1x,i3,1x,a10,1x,i2,x,a10,1x,i2)') & write(6,'(a10,1x,f8.4,1x,a10,1x,i4,1x,a10,1x,i3,1x,a10,1x,i2,x,a10,1x,i2)') &
'theTime',theTime,'theInc',theInc,'theCycle',theCycle,'theLovl',theLovl,& 'theTime',theTime,'theInc',theInc,'theCycle',theCycle,'theLovl',theLovl,&
'mode',CPFEM_mode 'mode',CPFEM_mode
write(6,*) '#####################################'
endif endif
!
select case (CPFEM_mode) select case (CPFEM_mode)
case (1,2) ! regular computation (with aging of results if mode == 1) case (1,2) ! regular computation (with aging of results if mode == 1)
if (CPFEM_mode == 1) then ! age results at start of new increment if (CPFEM_mode == 1) then ! age results at start of new increment
CPFEM_Lp_old = CPFEM_Lp_new crystallite_F0 = crystallite_partionedF ! crystallite deformation (_subF is perturbed...)
CPFEM_Fp_old = CPFEM_Fp_new crystallite_Fp0 = crystallite_Fp ! crystallite plastic deformation
crystallite_Lp0 = crystallite_Lp ! crystallite plastic velocity
forall (i = 1:homogenization_maxNgrains,& forall (i = 1:homogenization_maxNgrains,&
j = 1:mesh_maxNips, & j = 1:mesh_maxNips, &
k = 1:mesh_NcpElems) & k = 1:mesh_NcpElems) &
constitutive_state_old(i,j,k)%p = constitutive_state_new(i,j,k)%p constitutive_state0(i,j,k)%p = constitutive_state(i,j,k)%p ! microstructure of crystallites
write (6,*) 'results aged.' write(6,'(a10,/,4(3(f10.3,x),/))') 'aged state',constitutive_state(1,1,1)%p/1e6
do j = 1,mesh_maxNips
do k = 1,mesh_NcpElems
if (homogenization_sizeState(j,k) > 0_pInt) &
homogenization_state0(j,k)%p = homogenization_state(j,k)%p ! internal state of homogenization scheme
enddo
enddo
endif endif
if (outdatedFFN1 .or. any(abs(ffn1 - CPFEM_ffn1_bar(:,:,CPFEM_in,cp_en)) > relevantStrain)& if (outdatedFFN1 .or. any(abs(ffn1 - materialpoint_F(:,:,CPFEM_in,cp_en)) > relevantStrain)) then
.and. parallelExecution) then if (.not. outdatedFFN1) write(6,'(a11,x,i5,x,i2,x,a10,/,3(3(f10.3,x),/))') 'outdated at',cp_en,CPFEM_in,'FFN1 now:',ffn1(:,1),ffn1(:,2),ffn1(:,3)
if (.not. outdatedFFN1) write(6,'(i5,x,i2,x,a10,/,3(3(f10.3,x),/))') cp_en,CPFEM_in,'FFN1 now:',ffn1(:,1),ffn1(:,2),ffn1(:,3)
outdatedFFN1 = .true. outdatedFFN1 = .true.
CPFEM_stress_bar(1:CPFEM_ngens,CPFEM_in,cp_en) = CPFEM_odd_stress CPFEM_cs(1:CPFEM_ngens,CPFEM_in,cp_en) = CPFEM_odd_stress
CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = CPFEM_odd_jacobian*math_identity2nd(CPFEM_ngens) CPFEM_dcsde(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = CPFEM_odd_jacobian*math_identity2nd(CPFEM_ngens)
else else
if (.not. parallelExecution) then if (.not. parallelExecution) then
CPFEM_execution_elem(1) = cp_en FEsolving_execElem(1) = cp_en
CPFEM_execution_elem(2) = cp_en FEsolving_execElem(2) = cp_en
CPFEM_execution_IP(1,cp_en) = CPFEM_in FEsolving_execIP(1,cp_en) = CPFEM_in
CPFEM_execution_IP(2,cp_en) = CPFEM_in FEsolving_execIP(2,cp_en) = CPFEM_in
CPFEM_Temperature(CPFEM_in,cp_en) = Temperature call materialpoint_stressAndItsTangent(CPFEM_updateJaco, CPFEM_dt)
CPFEM_ffn_bar(:,:,CPFEM_in,cp_en) = ffn call materialpoint_postResults(CPFEM_dt)
CPFEM_ffn1_bar(:,:,CPFEM_in,cp_en) = ffn1
call CPFEM_MaterialPoint(CPFEM_updateJaco, CPFEM_dt)
elseif (.not. CPFEM_calc_done) then elseif (.not. CPFEM_calc_done) then
call CPFEM_MaterialPoint(CPFEM_updateJaco, CPFEM_dt) ! parallel execution inside call materialpoint_stressAndItsTangent(CPFEM_updateJaco, CPFEM_dt) ! parallel execution inside
call materialpoint_postResults(CPFEM_dt)
CPFEM_calc_done = .true. CPFEM_calc_done = .true.
endif endif
! translate from P and dP/dF to CS and dCS/dE ! translate from P and dP/dF to CS and dCS/dE
Kirchhoff_bar = math_mul33x33(CPFEM_PK1_bar(:,:,CPFEM_in, cp_en),transpose(CPFEM_ffn1_bar(:,:,CPFEM_in, cp_en))) Kirchhoff = math_mul33x33(materialpoint_P(:,:,CPFEM_in, cp_en),transpose(materialpoint_F(:,:,CPFEM_in, cp_en)))
J_inverse = 1.0_pReal/math_det3x3(CPFEM_ffn1_bar(:,:,CPFEM_in, cp_en)) J_inverse = 1.0_pReal/math_det3x3(materialpoint_F(:,:,CPFEM_in, cp_en))
CPFEM_stress_bar(1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel33to6(J_inverse*Kirchhoff_bar) CPFEM_cs(1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel33to6(J_inverse*Kirchhoff)
!
H_bar = 0.0_pReal H = 0.0_pReal
forall(i=1:3,j=1:3,k=1:3,l=1:3,m=1:3,n=1:3) & forall(i=1:3,j=1:3,k=1:3,l=1:3,m=1:3,n=1:3) &
H_bar(i,j,k,l) = H_bar(i,j,k,l) + & H(i,j,k,l) = H(i,j,k,l) + &
CPFEM_ffn1_bar(j,m,CPFEM_in,cp_en) * & materialpoint_F(j,m,CPFEM_in,cp_en) * &
CPFEM_ffn1_bar(l,n,CPFEM_in,cp_en) * & materialpoint_F(l,n,CPFEM_in,cp_en) * &
CPFEM_dPdF_bar(i,m,k,n,CPFEM_in,cp_en) - & materialpoint_dPdF(i,m,k,n,CPFEM_in,cp_en) - &
math_I3(j,l)*CPFEM_ffn1_bar(i,m,CPFEM_in,cp_en)*CPFEM_PK1_bar(k,m,CPFEM_in,cp_en) + & math_I3(j,l)*materialpoint_F(i,m,CPFEM_in,cp_en)*materialpoint_P(k,m,CPFEM_in,cp_en) + &
0.5_pReal*(math_I3(i,k)*Kirchhoff_bar(j,l) + math_I3(j,l)*Kirchhoff_bar(i,k) + & 0.5_pReal*(math_I3(i,k)*Kirchhoff(j,l) + math_I3(j,l)*Kirchhoff(i,k) + &
math_I3(i,l)*Kirchhoff_bar(j,k) + math_I3(j,k)*Kirchhoff_bar(i,l)) math_I3(i,l)*Kirchhoff(j,k) + math_I3(j,k)*Kirchhoff(i,l))
forall(i=1:3,j=1:3,k=1:3,l=1:3) & forall(i=1:3,j=1:3,k=1:3,l=1:3) &
H_bar_sym(i,j,k,l)= 0.25_pReal*(H_bar(i,j,k,l) +H_bar(j,i,k,l) +H_bar(i,j,l,k) +H_bar(j,i,l,k)) 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_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel3333to66(J_inverse*H_bar) CPFEM_dcsde(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel3333to66(J_inverse*H)
endif endif
case (3) ! collect and return odd result case (3) ! collect and return odd result
CPFEM_Temperature(CPFEM_in,cp_en) = Temperature materialpoint_Temperature(CPFEM_in,cp_en) = Temperature
CPFEM_ffn_bar(:,:,CPFEM_in,cp_en) = ffn materialpoint_F0(:,:,CPFEM_in,cp_en) = ffn
CPFEM_ffn1_bar(:,:,CPFEM_in,cp_en) = ffn1 materialpoint_F(:,:,CPFEM_in,cp_en) = ffn1
CPFEM_stress_bar(1:CPFEM_ngens,CPFEM_in,cp_en) = CPFEM_odd_stress CPFEM_cs(1:CPFEM_ngens,CPFEM_in,cp_en) = CPFEM_odd_stress
CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = CPFEM_odd_jacobian*math_identity2nd(CPFEM_ngens) CPFEM_dcsde(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = CPFEM_odd_jacobian*math_identity2nd(CPFEM_ngens)
CPFEM_calc_done = .false. CPFEM_calc_done = .false.
case (4) ! do nothing since we can recycle the former results (MARC specialty) case (4) ! do nothing since we can recycle the former results (MARC specialty)
case (5) ! record consistent tangent at beginning of new increment (while recycling) case (5) ! record consistent tangent at beginning of new FE increment (while recycling)
CPFEM_jaco_knownGood = CPFEM_jaco_bar CPFEM_dcsde_knownGood = CPFEM_dcsde
case (6) ! restore consistent tangent after cutback case (6) ! restore consistent tangent after FE cutback
CPFEM_jaco_bar = CPFEM_jaco_knownGood CPFEM_dcsde = CPFEM_dcsde_knownGood
end select end select
!
! return the local stress and the jacobian from storage ! return the local stress and the jacobian from storage
CPFEM_stress(1:CPFEM_ngens) = CPFEM_stress_bar(1:CPFEM_ngens,CPFEM_in,cp_en) CPFEM_stress(1:CPFEM_ngens) = CPFEM_cs(1:CPFEM_ngens,CPFEM_in,cp_en)
CPFEM_jaco(1:CPFEM_ngens,1:CPFEM_ngens) = CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) CPFEM_jaco(1:CPFEM_ngens,1:CPFEM_ngens) = CPFEM_dcsdE(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en)
!
return
!
END SUBROUTINE
!
!
!**********************************************************
!*** calculate the material point behaviour ***
!**********************************************************
SUBROUTINE CPFEM_MaterialPoint(&
updateJaco,& ! flag to initiate Jacobian updating
CPFEM_dt) ! Time increment (dt)
!
use prec
use debug
use math, only: math_pDecomposition,math_RtoEuler,inDeg
use IO, only: IO_error
use mesh, only: mesh_element, mesh_NcpElems, FE_Nips
use material, only: homogenization_Ngrains,material_phase,material_volfrac
use constitutive
implicit none
!
logical, intent(in) :: updateJaco
real(pReal), intent(in) :: CPFEM_dt
integer(pInt) g,i,e
logical error
real(pReal), dimension(3,3) :: U,R
!$OMP PARALLEL DO
do e = CPFEM_execution_elem(1),CPFEM_execution_elem(2) ! iterate over elements to be processed
do i = CPFEM_execution_IP(1,e),CPFEM_execution_IP(2,e) ! iterate over IPs of this element to be processed
forall (g = 1:homogenization_Ngrains(mesh_element(3,e))) ! number of grains of this homogenization
CPFEM_ffn(:,:,g,i,e) = CPFEM_ffn_bar(:,:,i,e) ! Taylor homogenization (why not using former ffn1??)
CPFEM_ffn1(:,:,g,i,e) = CPFEM_ffn1_bar(:,:,i,e) ! Taylor homogenization
end forall
enddo
enddo
!$OMP END PARALLEL DO
call SingleCrystallite(updateJaco,CPFEM_dt)
!******************************************************************************************************
! check convergence of homogenization if needed
!******************************************************************************************************
! calculate average quantities per ip and post results
!$OMP PARALLEL DO
do e = CPFEM_execution_elem(1),CPFEM_execution_elem(2) ! iterate over elements to be processed
do i = CPFEM_execution_IP(1,e),CPFEM_execution_IP(2,e) ! iterate over IPs of this element to be processed
CPFEM_PK1_bar(:,:,i,e) = sum(CPFEM_PK1(:,:,:,i,e),3)/homogenization_Ngrains(mesh_element(3,e))
if (updateJaco) &
CPFEM_dPdF_bar(:,:,:,:,i,e) = &
sum(CPFEM_dPdF(:,:,:,:,:,i,e),5)/homogenization_Ngrains(mesh_element(3,e)) ! add up crystallite stiffnesses (may have "holes" corresponding to former avg tangent)
do g = 1,homogenization_Ngrains(mesh_element(3,e))
call math_pDecomposition(CPFEM_Fe_new(:,:,g,i,e),U,R,error) ! polar decomposition
if (error) call IO_error(650,e,i,g)
CPFEM_results(1,g,i,e) = material_phase(g,i,e)
CPFEM_results(2,g,i,e) = material_volFrac(g,i,e)
CPFEM_results(3:5,g,i,e) = math_RtoEuler(transpose(R))*inDeg ! orientation
enddo
enddo
enddo
!$OMP END PARALLEL DO
return
END SUBROUTINE
!********************************************************************
! Calculates the stress and jacobi (if wanted) for all or a single component
!********************************************************************
subroutine SingleCrystallite(&
updateJaco,& ! update of Jacobian required
dt) ! time increment
use prec, only: pReal,pInt,pert_Fg,subStepMin, nCutback
use debug
use math
use IO, only: IO_error
use mesh, only: mesh_element, FE_Nips
use material, only: homogenization_Ngrains
use constitutive
implicit none
character (len=128) msg
logical updateJaco, allConverged
real(preal) dt
real(pReal), dimension(3,3) :: Fg_pert,Lp_pert, P_pert, Fp_pert, Fe_pert
real(pReal), dimension(6) :: Tstar_v
real(pReal), dimension(constitutive_maxSizeState) :: state
integer(pInt) g,i,e,k,l,iOuter,mySizeState
!$OMP PARALLEL DO
do e = CPFEM_execution_elem(1),CPFEM_execution_elem(2) ! iterate over elements to be processed
do i = CPFEM_execution_IP(1,e),CPFEM_execution_IP(2,e) ! iterate over IPs of this element to be processed
forall (g = 1:homogenization_Ngrains(mesh_element(3,e))) ! number of grains of this homogenization
crystallite_converged(g,i,e) = .false.
constitutive_state_new(g,i,e)%p = constitutive_state_old(g,i,e)%p
CPFEM_Lp_new(:,:,g,i,e) = CPFEM_Lp_old(:,:,g,i,e)
end forall
end do
end do
!$OMP END PARALLEL DO
iOuter = 0_pInt
allConverged = .false.
do while (.not. allConverged)
iOuter = iOuter + 1_pInt ! count state integation loops
if (iOuter > nOuter) call IO_error(600) ! too many loops required --> croak
!$OMP PARALLEL DO
do e = CPFEM_execution_elem(1),CPFEM_execution_elem(2) ! iterate over elements to be processed
do i = CPFEM_execution_IP(1,e),CPFEM_execution_IP(2,e) ! iterate over IPs of this element to be processed
do g = 1,homogenization_Ngrains(mesh_element(3,e)) ! number of grains of this homogenization
if (.not. crystallite_converged(g,i,e)) then
call integrateStress(msg,CPFEM_Tstar_v(:,g,i,e),CPFEM_PK1(:,:,g,i,e), &
CPFEM_Fp_new(:,:,g,i,e),CPFEM_Fe_new(:,:,g,i,e),CPFEM_Lp_new(:,:,g,i,e), &
CPFEM_ffn1(:,:,g,i,e),dt,g,i,e)
if (msg /= 'ok') call IO_error(610,e,i,g,msg)
endif
end do
end do
end do
!$OMP END PARALLEL DO
allConverged = .true. ! assume best case
!$OMP PARALLEL DO
do e = CPFEM_execution_elem(1),CPFEM_execution_elem(2) ! iterate over elements to be processed
do i = CPFEM_execution_IP(1,e),CPFEM_execution_IP(2,e) ! iterate over IPs of this element to be processed
do g = 1,homogenization_Ngrains(mesh_element(3,e)) ! number of grains of this homogenization
if (crystallite_converged(g,i,e)) cycle ! this one is already fine
if (integrateState(CPFEM_Tstar_v(:,g,i,e),dt,g,i,e)) then ! state integration now converged?
crystallite_converged(g,i,e) = .true.
!$OMP CRITICAL (out)
debug_OuterLoopDistribution(iOuter) = debug_OuterLoopDistribution(iOuter)+1
!$OMP END CRITICAL (out)
else
allConverged = .false. ! this one requires additional round...
endif
end do
end do
end do
!$OMP END PARALLEL DO
end do ! all crystallites converged
!$OMP PARALLEL DO
do e = CPFEM_execution_elem(1),CPFEM_execution_elem(2) ! iterate over elements to be processed
do i = CPFEM_execution_IP(1,e),CPFEM_execution_IP(2,e) ! iterate over IPs of this element to be processed
forall (g = 1:homogenization_Ngrains(mesh_element(3,e))) & ! number of grains of this homogenization
CPFEM_results(CPFEM_Nresults+1:CPFEM_Nresults+constitutive_sizePostResults(g,i,e),g,i,e) = &
constitutive_postResults(CPFEM_Tstar_v(:,g,i,e),CPFEM_Temperature(i,e),dt,g,i,e)
end do
end do
!$OMP END PARALLEL DO
if(updateJaco) then ! Jacobian required
!$OMP CRITICAL (write2out)
if (debugger) write (6,*) 'Jacobian calc'
!$OMP END CRITICAL (write2out)
!$OMP PARALLEL DO
do e = CPFEM_execution_elem(1),CPFEM_execution_elem(2) ! iterate over elements to be processed
do i = CPFEM_execution_IP(1,e),CPFEM_execution_IP(2,e) ! iterate over IPs of this element to be processed
do g = 1,homogenization_Ngrains(mesh_element(3,e)) ! number of grains of this homogenization
mySizeState = constitutive_sizeState(g,i,e) ! number of state variables for this grain
state(1:mySizeState) = constitutive_state_new(g,i,e)%p ! remember unperturbed, converged state
do k = 1,3 ! perturbation...
do l = 1,3 ! ...components
Fg_pert = CPFEM_ffn1(:,:,g,i,e) ! initialize perturbed Fg
Fg_pert(k,l) = Fg_pert(k,l) + pert_Fg ! perturb single component
Lp_pert = CPFEM_Lp_new(:,:,g,i,e) ! initialize Lp
Fp_pert = CPFEM_Fp_new(:,:,g,i,e) ! initialize Fp
constitutive_state_new(g,i,e)%p = state(1:mySizeState) ! initial guess from end of time step
crystallite_converged(g,i,e) = .false.
iOuter = 0_pInt
do while(.not. crystallite_converged(g,i,e) .and. iOuter < nOuter)
iOuter = iOuter + 1_pInt
call integrateStress(msg,Tstar_v,P_pert,Fp_pert,Fe_pert,Lp_pert, Fg_pert,dt,g,i,e)
if (msg /= 'ok') exit
crystallite_converged(g,i,e) = integrateState(Tstar_v,dt,g,i,e)
end do
if (crystallite_converged(g,i,e)) &
CPFEM_dPdF(:,:,k,l,g,i,e) = (P_pert-CPFEM_PK1(:,:,g,i,e))/pert_Fg ! constructing tangent dP_ij/dFg_kl only if valid forward difference
!$OMP CRITICAL (out)
debug_OuterLoopDistribution(iOuter) = debug_OuterLoopDistribution(iOuter)+1
!$OMP END CRITICAL (out)
end do
end do
constitutive_state_new(g,i,e)%p = state(1:mySizeState) ! restore solution
end do
end do
end do
!$OMP END PARALLEL DO
endif
return return
end subroutine end subroutine
!********************************************************************
! Update the state for a single component
!********************************************************************
function integrateState(&
Tstar_v,& ! stress
dt,& ! time increment
g,& ! grain number
i,& ! integration point number
e& ! element number
)
use prec, only: pReal,pInt,pLongInt,reltol_Outer
use constitutive, only: constitutive_dotState,constitutive_sizeDotState,&
constitutive_state_old,constitutive_state_new
use debug
logical integrateState
integer(pLongInt) tick,tock,tickrate,maxticks
integer(pInt) g,i,e,mySize
real(pReal), dimension(6) :: Tstar_v
real(pReal) dt
real(pReal), dimension(constitutive_sizeDotState(g,i,e)) :: residuum
mySize = constitutive_sizeDotState(g,i,e)
call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
residuum = constitutive_state_new(g,i,e)%p(1:mySize) - constitutive_state_old(g,i,e)%p(1:mySize) - &
dt*constitutive_dotState(Tstar_v,CPFEM_Temperature(i,e),g,i,e) ! residuum from evolution of microstructure
call system_clock(count=tock,count_rate=tickrate,count_max=maxticks)
debug_cumDotStateCalls = debug_cumDotStateCalls + 1_pInt
debug_cumDotStateTicks = debug_cumDotStateTicks + tock-tick
if (tock < tick) debug_cumDotStateTicks = debug_cumDotStateTicks + maxticks
constitutive_state_new(g,i,e)%p(1:mySize) = constitutive_state_new(g,i,e)%p(1:mySize) - residuum ! update of microstructure
integrateState = maxval(abs(residuum/constitutive_state_new(g,i,e)%p(1:mySize)),&
constitutive_state_new(g,i,e)%p(1:mySize) /= 0.0_pReal) < reltol_Outer
return
end function
!********************************************************************
! Calculates the stress for a single component
!********************************************************************
!***********************************************************************
!*** calculation of stress (P), stiffness (dPdF), ***
!*** and announcement of any ***
!*** acceleration of the Newton-Raphson correction ***
!***********************************************************************
subroutine integrateStress(&
msg,& ! return message
Tstar_v,& ! Stress vector
P,& ! first PK stress
Fp_new,& ! new plastic deformation gradient
Fe_new,& ! new "elastic" deformation gradient
Lp,& ! plastic velocity gradient
!
Fg_new,& ! new global deformation gradient
dt,& ! time increment
g,& ! grain number
i,& ! integration point number
e) ! element number
use prec, only: pReal,pInt,pert_Fg,subStepMin, nCutback
use debug
use constitutive, only: constitutive_state_new
use math
! use CPFEM
!
implicit none
!
character(len=*) msg
logical error,success
integer(pInt) e,i,g, nCutbacks, maxCutbacks
real(pReal) Temperature
real(pReal) dt,dt_aim,subFrac,subStep,det
real(pReal), dimension(3,3) :: Lp,Lp_interpolated,inv
real(pReal), dimension(3,3) :: Fg_current,Fg_new,Fg_aim,deltaFg
real(pReal), dimension(3,3) :: Fp_current,Fp_new
real(pReal), dimension(3,3) :: Fe_current,Fe_new
real(pReal), dimension(3,3) :: P
real(pReal), dimension(6) :: Tstar_v
deltaFg = Fg_new - CPFEM_ffn(:,:,g,i,e)
subFrac = 0.0_pReal
subStep = 1.0_pReal
nCutbacks = 0_pInt
maxCutbacks = 0_pInt
Fg_current = CPFEM_ffn(:,:,g,i,e) ! initialize to start of inc
Fp_current = CPFEM_Fp_old(:,:,g,i,e)
call math_invert3x3(Fp_current,inv,det,error)
Fe_current = math_mul33x33(Fg_current,inv)
success = .false. ! pretend cutback
dt_aim = 0.0_pReal ! prevent initial Lp interpolation
Temperature = CPFEM_Temperature(i,e)
! begin the cutback loop
do while (subStep > subStepMin) ! continue until finished or too much cut backing
if (success) then ! wind forward
Fg_current = Fg_aim
Fe_current = Fe_new
Fp_current = Fp_new
elseif (dt_aim > 0.0_pReal) then
call math_invert3x3(Fg_aim,inv,det,error) ! inv of Fg_aim
Lp_interpolated = 0.5_pReal*Lp + &
0.5_pReal*(math_I3 - math_mul33x33(Fp_current,&
math_mul33x33(inv,Fe_current)))/dt_aim ! interpolate Lp and L
if (debugger) then
!$OMP CRITICAL (write2out)
write (6,*) 'Lp interpolation'
write (6,'(a,/,3(3(f12.7,x)/))') 'from',Lp(1:3,:)
write (6,'(a,/,3(3(f12.7,x)/))') 'to',Lp_interpolated(1:3,:)
!$OMP END CRITICAL (write2out)
endif
Lp = Lp_interpolated
endif
!
Fg_aim = Fg_current + subStep*deltaFg ! aim for Fg
dt_aim = subStep*dt ! aim for dt
if (debugger) then
!$OMP CRITICAL (write2out)
write (6,*) 'using these values'
write (6,'(a,/,3(4(f9.3,x)/))') 'state new / MPa',constitutive_state_new(g,i,e)%p/1e6_pReal
write (6,'(a,/,3(3(f12.7,x)/))') 'Fe current',Fe_current(1:3,:)
write (6,'(a,/,3(3(f12.7,x)/))') 'Fp current',Fp_current(1:3,:)
write (6,'(a,/,3(3(f12.7,x)/))') 'Lp (old=new guess)',Lp(1:3,:)
write (6,'(a20,f,x,a2,x,f)') 'integrating from ',subFrac,'to',(subFrac+subStep)
!$OMP END CRITICAL (write2out)
endif
call TimeIntegration(msg,Lp,Fp_new,Fe_new,Tstar_v,P, Fg_aim,Fp_current,Temperature,dt_aim,g,i,e)
if (msg == 'ok') then
subFrac = subFrac + subStep
subStep = min(1.0_pReal-subFrac, subStep*2.0_pReal) ! accelerate
nCutbacks = 0_pInt ! reset cutback counter
success = .true. ! keep current Lp
else
nCutbacks = nCutbacks + 1 ! record additional cutback
maxCutbacks = max(nCutbacks,maxCutbacks) ! remember maximum number of cutbacks
subStep = subStep / 2.0_pReal ! cut time step in half
success = .false. ! force Lp interpolation
endif
enddo ! potential substepping
!
!$OMP CRITICAL (cutback)
debug_cutbackDistribution(min(nCutback,maxCutbacks)+1) = debug_cutbackDistribution(min(nCutback,maxCutbacks)+1)+1
!$OMP END CRITICAL (cutback)
return
end subroutine
!
!***********************************************************************
!*** fully-implicit two-level time integration ***
!*** based on a residuum in Lp and intermediate ***
!*** acceleration of the Newton-Raphson correction ***
!***********************************************************************
SUBROUTINE TimeIntegration(&
msg,& ! return message
Lpguess,& ! guess of plastic velocity gradient
Fp_new,& ! new plastic deformation gradient
Fe_new,& ! new "elastic" deformation gradient
Tstar_v,& ! Stress vector
P,& ! 1st PK stress (taken as initial guess if /= 0)
Fg_new,& ! new total def gradient
Fp_old,& ! former plastic def gradient
Temperature,& ! temperature
dt,& ! time increment
grain,& ! grain number
ip,& ! integration point number
cp_en & ! element number
)
use prec
use debug
use mesh, only: mesh_element
use constitutive, only: constitutive_microstructure,constitutive_homogenizedC,constitutive_LpAndItsTangent,&
constitutive_state_new
use math
use IO
implicit none
!
character(len=*) msg
logical failed
integer(pInt) cp_en, ip, grain
integer(pInt) iInner,dummy, i,j,k,l,m,n
integer(pLongInt) tick,tock,tickrate,maxticks
real(pReal) dt, Temperature, det, p_hydro, leapfrog,maxleap
real(pReal), dimension(6) :: Tstar_v
real(pReal), dimension(9,9) :: dLp,dTdLp,dRdLp,invdRdLp,eye2
real(pReal), dimension(6,6) :: C_66
real(pReal), dimension(3,3) :: Fg_new,Fp_new,invFp_new,Fp_old,invFp_old,Fe_new
real(pReal), dimension(3,3) :: P
real(pReal), dimension(3,3) :: Lp,Lpguess,Lpguess_old,Rinner,Rinner_old,A,B,BT,AB,BTA
real(pReal), dimension(3,3,3,3) :: C
msg = 'ok' ! error-free so far
eye2 = math_identity2nd(9)
call math_invert3x3(Fp_old,invFp_old,det,failed) ! inversion of Fp_old
if (failed) then
msg = 'inversion Fp_old'
return
endif
A = math_mul33x33(transpose(invFp_old), math_mul33x33(transpose(Fg_new),math_mul33x33(Fg_new,invFp_old)))
!$OMP CRITICAL (write2out)
if (debugger) write (6,'(a,/,3(3(f12.7,x)/))') 'Fg to be calculated',Fg_new
!$OMP END CRITICAL (write2out)
call constitutive_microstructure(Temperature,grain,ip,cp_en)
C_66 = constitutive_homogenizedC(grain,ip,cp_en)
C = math_Mandel66to3333(C_66) ! 4th rank elasticity tensor
iInner = 0_pInt
leapfrog = 1.0_pReal ! correction as suggested by invdRdLp-step
maxleap = 1024.0_pReal ! preassign maximum acceleration level
Lpguess_old = Lpguess ! consider present Lpguess good
Inner: do ! inner iteration: Lp
iInner = iInner+1
if (iInner > nInner) then ! too many loops required
Lpguess = Lpguess_old ! do not trust the last update but resort to former one
msg = 'limit Inner iteration'
return
endif
B = math_i3 - dt*Lpguess
BT = transpose(B)
AB = math_mul33x33(A,B)
BTA = math_mul33x33(BT,A)
Tstar_v = 0.5_pReal*math_mul66x6(C_66,math_mandel33to6(math_mul33x33(BT,AB)-math_I3))
p_hydro=(Tstar_v(1)+Tstar_v(2)+Tstar_v(3))/3.0_pReal
forall(i=1:3) Tstar_v(i) = Tstar_v(i)-p_hydro ! subtract hydrostatic pressure
call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
call constitutive_LpAndItsTangent(Lp,dLp, Tstar_v,Temperature,grain,ip,cp_en)
call system_clock(count=tock,count_rate=tickrate,count_max=maxticks)
debug_cumLpCalls = debug_cumLpCalls + 1_pInt
debug_cumLpTicks = debug_cumLpTicks + tock-tick
if (tock < tick) debug_cumLpTicks = debug_cumLpTicks + maxticks
Rinner = Lpguess - Lp ! update current residuum
if (.not.(any(Rinner/=Rinner)) .and. & ! exclude any NaN in residuum
( ( maxval(abs(Rinner)) < abstol_Inner) .or. & ! below abs tol .or.
( any(abs(dt*Lpguess) > relevantStrain) .and. & ! worth checking? .and.
maxval(abs(Rinner/Lpguess),abs(dt*Lpguess) > relevantStrain) < reltol_Inner & ! below rel tol
) &
) &
) &
exit Inner ! convergence
!
! check for acceleration/deceleration in Newton--Raphson correction
!
if (any(Rinner/=Rinner) .and. & ! NaN occured at regular speed
leapfrog == 1.0) then
Lpguess = Lpguess_old ! restore known good guess
msg = 'NaN present' ! croak for cutback
return
elseif (leapfrog > 1.0_pReal .and. & ! at fast pace ?
(sum(Rinner*Rinner) > sum(Rinner_old*Rinner_old) .or. & ! worse residuum
sum(Rinner*Rinner_old) < 0.0_pReal) .or. & ! residuum changed sign (overshoot)
any(Rinner/=Rinner) ) then ! NaN
maxleap = 0.5_pReal * leapfrog ! limit next acceleration
leapfrog = 1.0_pReal ! grinding halt
else ! better residuum
dTdLp = 0.0_pReal ! calc dT/dLp
forall (i=1:3,j=1:3,k=1:3,l=1:3,m=1:3,n=1:3) &
dTdLp(3*(i-1)+j,3*(k-1)+l) = dTdLp(3*(i-1)+j,3*(k-1)+l) + &
C(i,j,l,n)*AB(k,n)+C(i,j,m,l)*BTA(m,k)
dTdLp = -0.5_pReal*dt*dTdLp
dRdLp = eye2 - math_mul99x99(dLp,dTdLp) ! calc dR/dLp
invdRdLp = 0.0_pReal
call math_invert(9,dRdLp,invdRdLp,dummy,failed) ! invert dR/dLp --> dLp/dR
if (failed) then
msg = 'inversion dR/dLp'
if (debugger) then
!$OMP CRITICAL (write2out)
write (6,*) msg
write (6,'(a,/,9(9(e9.3,x)/))') 'dRdLp', dRdLp(1:9,:)
write (6,'(a,/,3(4(f9.3,x)/))') 'state_new / MPa',constitutive_state_new(grain,ip,cp_en)%p/1e6_pReal
write (6,'(a,/,3(3(f12.7,x)/))') 'Lpguess',Lpguess(1:3,:)
write (6,'(a,/,3(3(e12.7,x)/))') 'Lp',Lp(1:3,:)
write (6,'(a,/,6(f9.3,x))') 'Tstar / MPa',Tstar_v/1e6_pReal
!$OMP END CRITICAL (write2out)
endif
return
endif
!
Rinner_old = Rinner ! remember current residuum
Lpguess_old = Lpguess ! remember current Lp guess
if (iInner > 1 .and. leapfrog < maxleap) leapfrog = 2.0_pReal * leapfrog ! accelerate if ok
endif
!
Lpguess = Lpguess_old ! start from current guess
Rinner = Rinner_old ! use current residuum
forall (i=1:3,j=1:3,k=1:3,l=1:3) & ! leapfrog to updated Lpguess
Lpguess(i,j) = Lpguess(i,j) - leapfrog*invdRdLp(3*(i-1)+j,3*(k-1)+l)*Rinner(k,l)
enddo Inner
!
!$OMP CRITICAL (in)
debug_InnerLoopDistribution(iInner) = debug_InnerLoopDistribution(iInner)+1
!$OMP END CRITICAL (in)
invFp_new = math_mul33x33(invFp_old,B)
call math_invert3x3(invFp_new,Fp_new,det,failed)
if (failed) then
msg = 'inversion Fp_new^-1'
return
endif
Fp_new = Fp_new*det**(1.0_pReal/3.0_pReal) ! regularize Fp by det = det(InvFp_new) !!
forall (i=1:3) Tstar_v(i) = Tstar_v(i) + p_hydro ! add hydrostatic component back
Fe_new = math_mul33x33(Fg_new,invFp_new) ! calc resulting Fe
P = math_mul33x33(Fe_new,math_mul33x33(math_Mandel6to33(Tstar_v),transpose(invFp_new))) ! first PK stress
return
!
END SUBROUTINE
!
END MODULE END MODULE
!############################################################## !##############################################################

13
trunk/FEsolving.f90
View File

@ -12,6 +12,8 @@
logical :: lastIncConverged = .false.,outdatedByNewInc = .false.,outdatedFFN1 = .false. logical :: lastIncConverged = .false.,outdatedByNewInc = .false.,outdatedFFN1 = .false.
logical :: symmetricSolver = .false. logical :: symmetricSolver = .false.
logical :: parallelExecution = .true. logical :: parallelExecution = .true.
integer(pInt), dimension(:,:), allocatable :: FEsolving_execIP
integer(pInt), dimension(2) :: FEsolving_execElem
CONTAINS CONTAINS
@ -26,7 +28,8 @@
implicit none implicit none
integer(pInt), parameter :: fileunit = 222 integer(pInt), parameter :: fileunit = 222
integer(pInt), dimension (1+2*2) :: pos integer(pInt), parameter :: maxNchunks = 2
integer(pInt), dimension(1+2*maxNchunks) :: positions
character(len=1024) line character(len=1024) line
if (IO_open_inputFile(fileunit)) then if (IO_open_inputFile(fileunit)) then
@ -34,11 +37,11 @@
rewind(fileunit) rewind(fileunit)
do do
read (fileunit,'(a1024)',END=100) line read (fileunit,'(a1024)',END=100) line
pos = IO_stringPos(line,1) positions = IO_stringPos(line,1)
if( IO_lc(IO_stringValue(line,pos,1)) == 'solver' ) then if( IO_lc(IO_stringValue(line,positions,1)) == 'solver' ) then
read (fileunit,'(a1024)',END=100) line ! Garbage line read (fileunit,'(a1024)',END=100) line ! Garbage line
pos = IO_stringPos(line,2) positions = IO_stringPos(line,2)
symmetricSolver = (IO_intValue(line,pos,2) /= 1_pInt) symmetricSolver = (IO_intValue(line,positions,2) /= 1_pInt)
exit exit
endif endif
enddo enddo

2
trunk/IO.f90
View File

@ -879,6 +879,8 @@ END FUNCTION
select case (ID) select case (ID)
case (650) case (650)
msg = 'Polar decomposition failed' msg = 'Polar decomposition failed'
case (600)
msg = 'Crystallite responds elastically'
case default case default
msg = 'Unknown warning number...' msg = 'Unknown warning number...'
end select end select

89
trunk/constitutive.f90
View File

@ -13,8 +13,10 @@ MODULE constitutive
use prec use prec
implicit none implicit none
type(p_vec), dimension(:,:,:), allocatable :: constitutive_state_old, & ! pointer array to old state variables of each grain type(p_vec), dimension(:,:,:), allocatable :: constitutive_state0, & ! pointer array to microstructure at start of FE inc
constitutive_state_new ! pointer array to new state variables of each grain 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)
integer(pInt), dimension(:,:,:), allocatable :: constitutive_sizeDotState, & ! size of dotState array integer(pInt), dimension(:,:,:), allocatable :: constitutive_sizeDotState, & ! size of dotState array
constitutive_sizeState, & ! size of state array per grain constitutive_sizeState, & ! size of state array per grain
constitutive_sizePostResults ! size of postResults array per grain constitutive_sizePostResults ! size of postResults array per grain
@ -44,7 +46,7 @@ subroutine constitutive_init()
use constitutive_dislobased use constitutive_dislobased
integer(pInt), parameter :: fileunit = 200 integer(pInt), parameter :: fileunit = 200
integer(pInt) e,i,g,myInstance integer(pInt) e,i,g,myInstance,myNgrains
if(.not. IO_open_file(fileunit,material_configFile)) call IO_error (100) ! corrupt config file if(.not. IO_open_file(fileunit,material_configFile)) call IO_error (100) ! corrupt config file
@ -54,44 +56,51 @@ subroutine constitutive_init()
close(fileunit) close(fileunit)
allocate(constitutive_state_old(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) allocate(constitutive_state0(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
allocate(constitutive_state_new(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) allocate(constitutive_partionedState0(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
allocate(constitutive_subState0(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
allocate(constitutive_state(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
allocate(constitutive_sizeDotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_sizeDotState = 0_pInt allocate(constitutive_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_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 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 do i = 1,FE_Nips(mesh_element(2,e)) ! loop over IPs
do g = 1,homogenization_Ngrains(mesh_element(3,e)) ! loop over grains do g = 1,myNgrains ! loop over grains
myInstance = phase_constitutionInstance(material_phase(g,i,e)) myInstance = phase_constitutionInstance(material_phase(g,i,e))
select case(phase_constitution(material_phase(g,i,e))) select case(phase_constitution(material_phase(g,i,e)))
case (constitutive_phenomenological_label) case (constitutive_phenomenological_label)
allocate(constitutive_state_old(g,i,e)%p(constitutive_phenomenological_sizeState(myInstance))) allocate(constitutive_state0(g,i,e)%p(constitutive_phenomenological_sizeState(myInstance)))
allocate(constitutive_state_new(g,i,e)%p(constitutive_phenomenological_sizeState(myInstance))) allocate(constitutive_partionedState0(g,i,e)%p(constitutive_phenomenological_sizeState(myInstance)))
constitutive_state_new(g,i,e)%p = constitutive_phenomenological_stateInit(myInstance) allocate(constitutive_subState0(g,i,e)%p(constitutive_phenomenological_sizeState(myInstance)))
constitutive_state_old(g,i,e)%p = constitutive_phenomenological_stateInit(myInstance) allocate(constitutive_state(g,i,e)%p(constitutive_phenomenological_sizeState(myInstance)))
constitutive_state0(g,i,e)%p = constitutive_phenomenological_stateInit(myInstance)
constitutive_sizeDotState(g,i,e) = constitutive_phenomenological_sizeDotState(myInstance) constitutive_sizeDotState(g,i,e) = constitutive_phenomenological_sizeDotState(myInstance)
constitutive_sizeState(g,i,e) = constitutive_phenomenological_sizeState(myInstance) constitutive_sizeState(g,i,e) = constitutive_phenomenological_sizeState(myInstance)
constitutive_sizePostResults(g,i,e) = constitutive_phenomenological_sizePostResults(myInstance) constitutive_sizePostResults(g,i,e) = constitutive_phenomenological_sizePostResults(myInstance)
case (constitutive_j2_label) case (constitutive_j2_label)
allocate(constitutive_state_old(g,i,e)%p(constitutive_j2_sizeState(myInstance))) allocate(constitutive_state0(g,i,e)%p(constitutive_j2_sizeState(myInstance)))
allocate(constitutive_state_new(g,i,e)%p(constitutive_j2_sizeState(myInstance))) allocate(constitutive_partionedState0(g,i,e)%p(constitutive_j2_sizeState(myInstance)))
constitutive_state_new(g,i,e)%p = constitutive_j2_stateInit(myInstance) allocate(constitutive_subState0(g,i,e)%p(constitutive_j2_sizeState(myInstance)))
constitutive_state_old(g,i,e)%p = constitutive_j2_stateInit(myInstance) allocate(constitutive_state(g,i,e)%p(constitutive_j2_sizeState(myInstance)))
constitutive_state0(g,i,e)%p = constitutive_j2_stateInit(myInstance)
constitutive_sizeDotState(g,i,e) = constitutive_j2_sizeDotState(myInstance) constitutive_sizeDotState(g,i,e) = constitutive_j2_sizeDotState(myInstance)
constitutive_sizeState(g,i,e) = constitutive_j2_sizeState(myInstance) constitutive_sizeState(g,i,e) = constitutive_j2_sizeState(myInstance)
constitutive_sizePostResults(g,i,e) = constitutive_j2_sizePostResults(myInstance) constitutive_sizePostResults(g,i,e) = constitutive_j2_sizePostResults(myInstance)
case (constitutive_dislobased_label) case (constitutive_dislobased_label)
allocate(constitutive_state_old(g,i,e)%p(constitutive_dislobased_sizeState(myInstance))) allocate(constitutive_state0(g,i,e)%p(constitutive_dislobased_sizeState(myInstance)))
allocate(constitutive_state_new(g,i,e)%p(constitutive_dislobased_sizeState(myInstance))) allocate(constitutive_partionedState0(g,i,e)%p(constitutive_dislobased_sizeState(myInstance)))
constitutive_state_new(g,i,e)%p = constitutive_dislobased_stateInit(myInstance) allocate(constitutive_subState0(g,i,e)%p(constitutive_dislobased_sizeState(myInstance)))
constitutive_state_old(g,i,e)%p = constitutive_dislobased_stateInit(myInstance) allocate(constitutive_state(g,i,e)%p(constitutive_dislobased_sizeState(myInstance)))
constitutive_state0(g,i,e)%p = constitutive_dislobased_stateInit(myInstance)
constitutive_sizeDotState(g,i,e) = constitutive_dislobased_sizeDotState(myInstance) constitutive_sizeDotState(g,i,e) = constitutive_dislobased_sizeDotState(myInstance)
constitutive_sizeState(g,i,e) = constitutive_dislobased_sizeState(myInstance) constitutive_sizeState(g,i,e) = constitutive_dislobased_sizeState(myInstance)
constitutive_sizePostResults(g,i,e) = constitutive_dislobased_sizePostResults(myInstance) constitutive_sizePostResults(g,i,e) = constitutive_dislobased_sizePostResults(myInstance)
case default case default
call IO_error(200,material_phase(g,i,e)) ! unknown constitution call IO_error(200,material_phase(g,i,e)) ! unknown constitution
end select end select
constitutive_partionedState0(g,i,e)%p = constitutive_state0(g,i,e)%p
enddo enddo
enddo enddo
enddo enddo
@ -99,6 +108,20 @@ subroutine constitutive_init()
constitutive_maxSizeState = maxval(constitutive_sizeState) constitutive_maxSizeState = maxval(constitutive_sizeState)
constitutive_maxSizePostResults = maxval(constitutive_sizePostResults) constitutive_maxSizePostResults = maxval(constitutive_sizePostResults)
write(6,*)
write(6,*) '<<<+- constitutive init -+>>>'
write(6,*)
write(6,'(a32,x,7(i5,x))') 'constitutive_state0: ', shape(constitutive_state0)
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_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)
write(6,*)
write(6,'(a32,x,7(i5,x))') 'maxSizeState: ', constitutive_maxSizeState
write(6,'(a32,x,7(i5,x))') 'maxSizePostResults: ', constitutive_maxSizePostResults
return return
end subroutine end subroutine
@ -126,11 +149,11 @@ function constitutive_homogenizedC(ipc,ip,el)
select case (phase_constitution(material_phase(ipc,ip,el))) select case (phase_constitution(material_phase(ipc,ip,el)))
case (constitutive_phenomenological_label) case (constitutive_phenomenological_label)
constitutive_homogenizedC = constitutive_phenomenological_homogenizedC(constitutive_state_new,ipc,ip,el) constitutive_homogenizedC = constitutive_phenomenological_homogenizedC(constitutive_state,ipc,ip,el)
case (constitutive_j2_label) case (constitutive_j2_label)
constitutive_homogenizedC = constitutive_j2_homogenizedC(constitutive_state_new,ipc,ip,el) constitutive_homogenizedC = constitutive_j2_homogenizedC(constitutive_state,ipc,ip,el)
case (constitutive_dislobased_label) case (constitutive_dislobased_label)
constitutive_homogenizedC = constitutive_dislobased_homogenizedC(constitutive_state_new,ipc,ip,el) constitutive_homogenizedC = constitutive_dislobased_homogenizedC(constitutive_state,ipc,ip,el)
end select end select
@ -161,11 +184,11 @@ real(pReal) Temperature
select case (phase_constitution(material_phase(ipc,ip,el))) select case (phase_constitution(material_phase(ipc,ip,el)))
case (constitutive_phenomenological_label) case (constitutive_phenomenological_label)
call constitutive_phenomenological_microstructure(Temperature,constitutive_state_new,ipc,ip,el) call constitutive_phenomenological_microstructure(Temperature,constitutive_state,ipc,ip,el)
case (constitutive_j2_label) case (constitutive_j2_label)
call constitutive_j2_microstructure(Temperature,constitutive_state_new,ipc,ip,el) call constitutive_j2_microstructure(Temperature,constitutive_state,ipc,ip,el)
case (constitutive_dislobased_label) case (constitutive_dislobased_label)
call constitutive_dislobased_microstructure(Temperature,constitutive_state_new,ipc,ip,el) call constitutive_dislobased_microstructure(Temperature,constitutive_state,ipc,ip,el)
end select end select
@ -201,11 +224,11 @@ subroutine constitutive_LpAndItsTangent(Lp,dLp_dTstar, Tstar_v,Temperature,ipc,i
select case (phase_constitution(material_phase(ipc,ip,el))) select case (phase_constitution(material_phase(ipc,ip,el)))
case (constitutive_phenomenological_label) case (constitutive_phenomenological_label)
call constitutive_phenomenological_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state_new,ipc,ip,el) call constitutive_phenomenological_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_j2_label) case (constitutive_j2_label)
call constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state_new,ipc,ip,el) call constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_dislobased_label) case (constitutive_dislobased_label)
call constitutive_dislobased_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state_new,ipc,ip,el) call constitutive_dislobased_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el)
end select end select
@ -241,11 +264,11 @@ function constitutive_dotState(Tstar_v,Temperature,ipc,ip,el)
select case (phase_constitution(material_phase(ipc,ip,el))) select case (phase_constitution(material_phase(ipc,ip,el)))
case (constitutive_phenomenological_label) case (constitutive_phenomenological_label)
constitutive_dotState = constitutive_phenomenological_dotState(Tstar_v,Temperature,constitutive_state_new,ipc,ip,el) constitutive_dotState = constitutive_phenomenological_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_j2_label) case (constitutive_j2_label)
constitutive_dotState = constitutive_j2_dotState(Tstar_v,Temperature,constitutive_state_new,ipc,ip,el) constitutive_dotState = constitutive_j2_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_dislobased_label) case (constitutive_dislobased_label)
constitutive_dotState = constitutive_dislobased_dotState(Tstar_v,Temperature,constitutive_state_new,ipc,ip,el) constitutive_dotState = constitutive_dislobased_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el)
end select end select
return return
@ -278,11 +301,11 @@ pure function constitutive_postResults(Tstar_v,Temperature,dt,ipc,ip,el)
constitutive_postResults = 0.0_pReal constitutive_postResults = 0.0_pReal
select case (phase_constitution(material_phase(ipc,ip,el))) select case (phase_constitution(material_phase(ipc,ip,el)))
case (constitutive_phenomenological_label) case (constitutive_phenomenological_label)
constitutive_postResults = constitutive_phenomenological_postResults(Tstar_v,Temperature,dt,constitutive_state_new,ipc,ip,el) constitutive_postResults = constitutive_phenomenological_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el)
case (constitutive_j2_label) case (constitutive_j2_label)
constitutive_postResults = constitutive_j2_postResults(Tstar_v,Temperature,dt,constitutive_state_new,ipc,ip,el) constitutive_postResults = constitutive_j2_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el)
case (constitutive_dislobased_label) case (constitutive_dislobased_label)
constitutive_postResults = constitutive_dislobased_postResults(Tstar_v,Temperature,dt,constitutive_state_new,ipc,ip,el) constitutive_postResults = constitutive_dislobased_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el)
end select end select

788
trunk/crystallite.f90 Normal file
View File

@ -0,0 +1,788 @@
!***************************************
!* Module: CRYSTALLITE *
!***************************************
!* contains: *
!* - _init *
!* - materialpoint_stressAndItsTangent *
!* - _partitionDeformation *
!* - _updateState *
!* - _averageStressAndItsTangent *
!* - _postResults *
!***************************************
MODULE crystallite
use prec, only: pReal,pInt
implicit none
!
! ****************************************************************
! *** General variables for the crystallite calculation ***
! ****************************************************************
integer(pInt), parameter :: crystallite_Nresults = 5_pInt ! phaseID, volfrac within this phase, Euler angles
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_Lp, & ! current plastic velocitiy grad (end of converged time step)
crystallite_F0, & ! def grad at start of FE inc
crystallite_Fp0, & ! plastic def grad at start of FE inc
crystallite_Lp0, & ! plastic velocitiy grad at start of FE inc
crystallite_partionedF, & ! def grad to be reached at end of homog inc
crystallite_partionedF0, & ! def grad at start of homog inc
crystallite_partionedFp0,& ! plastic def grad at start of homog inc
crystallite_partionedLp0,& ! plastic velocity grad at start of homog inc
crystallite_subF, & ! def grad to be reached at end of crystallite inc
crystallite_subF0, & ! def grad at start of crystallite inc
crystallite_subFp0,& ! plastic def grad at start of crystallite inc
crystallite_subLp0,& ! plastic velocity grad at start of crystallite inc
crystallite_P ! 1st Piola-Kirchhoff stress per grain
real(pReal), dimension (:,:,:,:), allocatable :: crystallite_Tstar_v ! 2nd Piola-Kirchhoff stress (vector) per grain
real(pReal), dimension (:,:,:,:,:,:,:),allocatable :: crystallite_dPdF, & ! individual dPdF per grain
crystallite_fallbackdPdF ! dPdF fallback for non-converged grains (elastic prediction)
real(pReal), dimension (:,:,:), allocatable :: crystallite_dt, & ! requested time increment of each grain
crystallite_subdt, & ! substepped time increment of each grain
crystallite_subFrac, & ! already calculated fraction of increment
crystallite_subStep, & ! size of next integration step
crystallite_Temperature ! Temp of each grain
logical, dimension (:,:,:), allocatable :: crystallite_localConstitution, & ! indicates this grain to have purely local constitutive law
crystallite_requested, & ! flag to request crystallite calculation
crystallite_onTrack, & ! flag to indicate ongoing calculation
crystallite_converged ! convergence flag
CONTAINS
!********************************************************************
! allocate and initialize per grain variables
!********************************************************************
subroutine crystallite_init()
use prec, only: pInt,pReal
use debug, only: debug_info,debug_reset
use math, only: math_I3,math_EulerToR
use FEsolving, only: FEsolving_execElem,FEsolving_execIP
use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips
use material, only: homogenization_Ngrains,homogenization_maxNgrains,&
material_EulerAngles,material_phase,phase_localConstitution
implicit none
integer(pInt) g,i,e, gMax,iMax,eMax, myNgrains
gMax = homogenization_maxNgrains
iMax = mesh_maxNips
eMax = mesh_NcpElems
allocate(crystallite_Fe(3,3,gMax,iMax,eMax)); crystallite_Fe = 0.0_pReal
allocate(crystallite_Fp(3,3,gMax,iMax,eMax)); crystallite_Fp = 0.0_pReal
allocate(crystallite_Lp(3,3,gMax,iMax,eMax)); crystallite_Lp = 0.0_pReal
allocate(crystallite_F0(3,3,gMax,iMax,eMax)); crystallite_F0 = 0.0_pReal
allocate(crystallite_Fp0(3,3,gMax,iMax,eMax)); crystallite_Fp0 = 0.0_pReal
allocate(crystallite_Lp0(3,3,gMax,iMax,eMax)); crystallite_Lp0 = 0.0_pReal
allocate(crystallite_partionedF(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_partionedLp0(3,3,gMax,iMax,eMax)); crystallite_partionedLp0 = 0.0_pReal
allocate(crystallite_subF(3,3,gMax,iMax,eMax)); crystallite_subF = 0.0_pReal
allocate(crystallite_subF0(3,3,gMax,iMax,eMax)); crystallite_subF0 = 0.0_pReal
allocate(crystallite_subFp0(3,3,gMax,iMax,eMax)); crystallite_subFp0 = 0.0_pReal
allocate(crystallite_subLp0(3,3,gMax,iMax,eMax)); crystallite_subLp0 = 0.0_pReal
allocate(crystallite_P(3,3,gMax,iMax,eMax)); crystallite_P = 0.0_pReal
allocate(crystallite_Tstar_v(6,gMax,iMax,eMax)); crystallite_Tstar_v = 0.0_pReal
allocate(crystallite_dPdF(3,3,3,3,gMax,iMax,eMax)); crystallite_dPdF = 0.0_pReal
allocate(crystallite_fallbackdPdF(3,3,3,3,gMax,iMax,eMax)); crystallite_fallbackdPdF = 0.0_pReal
allocate(crystallite_dt(gMax,iMax,eMax)); crystallite_dt = 0.0_pReal
allocate(crystallite_subdt(gMax,iMax,eMax)); crystallite_subdt = 0.0_pReal
allocate(crystallite_subFrac(gMax,iMax,eMax)); crystallite_subFrac = 0.0_pReal
allocate(crystallite_subStep(gMax,iMax,eMax)); crystallite_subStep = 0.0_pReal
allocate(crystallite_Temperature(gMax,iMax,eMax)); crystallite_Temperature = 0.0_pReal
allocate(crystallite_localConstitution(gMax,iMax,eMax));
allocate(crystallite_requested(gMax,iMax,eMax)); crystallite_requested = .false.
allocate(crystallite_onTrack(gMax,iMax,eMax)); crystallite_onTrack = .false.
allocate(crystallite_converged(gMax,iMax,eMax)); crystallite_converged = .true.
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over all cp elements
myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element
do g = 1,myNgrains
crystallite_Fp0(:,:,g,i,e) = math_EulerToR(material_EulerAngles(:,g,i,e)) ! plastic def gradient reflects init orientation
crystallite_F0(:,:,g,i,e) = math_I3
crystallite_partionedFp0(:,:,g,i,e) = crystallite_Fp0(:,:,g,i,e)
crystallite_partionedF0(:,:,g,i,e) = crystallite_F0(:,:,g,i,e)
crystallite_partionedF(:,:,g,i,e) = crystallite_F0(:,:,g,i,e)
crystallite_requested(g,i,e) = .true.
crystallite_localConstitution(g,i,e) = phase_localConstitution(material_phase(g,i,e))
enddo
enddo
enddo
!$OMP END PARALLEL DO
call crystallite_stressAndItsTangent(.true.) ! request elastic answers
crystallite_fallbackdPdF = crystallite_dPdF ! use initial elastic stiffness as fallback
! *** Output to MARC output file ***
!$OMP CRITICAL (write2out)
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_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_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_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_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_Temperature: ', shape(crystallite_Temperature)
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,*)
write(6,*) 'Number of non-local grains: ',count(.not. crystallite_localConstitution)
call flush(6)
!$OMP END CRITICAL (write2out)
call debug_info()
call debug_reset()
return
end subroutine
!********************************************************************
! calculate stress (P) and tangent (dPdF) for crystallites
!********************************************************************
subroutine crystallite_stressAndItsTangent(updateJaco)
use prec, only: pInt,pReal,subStepMin,nCryst
use debug
use IO, only: IO_warning
use math
use FEsolving, only: FEsolving_execElem, FEsolving_execIP
use mesh, only: mesh_element
use material, only: homogenization_Ngrains
use constitutive
implicit none
logical, intent(in) :: updateJaco
real(pReal), dimension(3,3) :: invFp,Fe_guess,PK2,myF,myFp,myFe,myLp,myP
real(pReal), dimension(constitutive_maxSizeState) :: myState
integer(pInt) crystallite_Niteration
integer(pInt) g,i,e,k,l, myNgrains, mySizeState
logical, dimension(2) :: doneAndHappy
! ------ initialize to starting condition ------
write (6,*)
write (6,*) 'Crystallite request from Materialpoint'
write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedF0 of 1 8 1',crystallite_partionedF0(1:3,:,1,8,1)
write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedFp0 of 1 8 1',crystallite_partionedFp0(1:3,:,1,8,1)
write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedF of 1 8 1',crystallite_partionedF(1:3,:,1,8,1)
write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedLp0 of 1 8 1',crystallite_partionedLp0(1:3,:,1,8,1)
!$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)) then ! initialize restoration point of ...
constitutive_subState0(g,i,e)%p = constitutive_partionedState0(g,i,e)%p ! ...microstructure
crystallite_subFp0(:,:,g,i,e) = crystallite_partionedFp0(:,:,g,i,e) ! ...plastic def grad
crystallite_subLp0(:,:,g,i,e) = crystallite_partionedLp0(:,:,g,i,e) ! ...plastic velocity grad
crystallite_subF0(:,:,g,i,e) = crystallite_partionedF0(:,:,g,i,e) ! ...def grad
crystallite_subFrac(g,i,e) = 0.0_pReal
crystallite_subStep(g,i,e) = 2.0_pReal
crystallite_onTrack(g,i,e) = .true.
crystallite_converged(g,i,e) = .false. ! pretend failed step of twice the required size
endif
enddo
enddo
enddo
!$OMP END PARALLEL DO
! ------ cutback loop ------
do while (any(crystallite_subStep(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) > subStepMin)) ! cutback loop for crystallites
write(6,*)
write(6,*) 'entering cutback at crystallite_stress'
if (any(.not. crystallite_converged(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) .and. & ! any non-converged grain
.not. crystallite_localConstitution(:,:,FEsolving_execELem(1):FEsolving_execElem(2))) ) & ! has non-local constitution?
crystallite_converged(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) = &
crystallite_converged(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) .and. &
crystallite_localConstitution(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) ! reset non-local grains' convergence status
!$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_converged(g,i,e)) then
crystallite_subFrac(g,i,e) = crystallite_subFrac(g,i,e) + crystallite_subStep(g,i,e)
crystallite_subStep(g,i,e) = min(1.0_pReal-crystallite_subFrac(g,i,e), 2.0_pReal * crystallite_subStep(g,i,e))
if (crystallite_subStep(g,i,e) > subStepMin) then
crystallite_subF0(:,:,g,i,e) = crystallite_subF(:,:,g,i,e) ! wind forward...
crystallite_subFp0(:,:,g,i,e) = crystallite_Fp(:,:,g,i,e) ! ...plastic def grad
crystallite_subLp0(:,:,g,i,e) = crystallite_Lp(:,:,g,i,e) ! ...plastic velocity gradient
constitutive_subState0(g,i,e)%p = constitutive_state(g,i,e)%p ! ...microstructure
endif
else
crystallite_subStep(g,i,e) = 0.5_pReal * crystallite_subStep(g,i,e) ! cut step in half and restore...
crystallite_Fp(:,:,g,i,e) = crystallite_subFp0(:,:,g,i,e) ! ...plastic def grad
crystallite_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
endif
crystallite_onTrack(g,i,e) = crystallite_subStep(g,i,e) > subStepMin ! still on track or already done (beyond repair)
if (crystallite_onTrack(g,i,e)) then ! specify task (according to substep)
crystallite_subF(:,:,g,i,e) = crystallite_subF0(:,:,g,i,e) + &
crystallite_subStep(g,i,e) * &
(crystallite_partionedF(:,:,g,i,e) - crystallite_partionedF0(:,:,g,i,e))
crystallite_subdt(g,i,e) = crystallite_subStep(g,i,e) * crystallite_dt(g,i,e)
crystallite_converged(g,i,e) = .false. ! start out non-converged
endif
enddo
enddo
enddo
!$OMP END PARALLEL DO
! ------ convergence loop for stress and state ------
crystallite_Niteration = 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)) &
) .and. crystallite_Niteration < nCryst) ! convergence loop for crystallite
crystallite_Niteration = crystallite_Niteration + 1
! --+>> stress integration <<+--
!
! incrementing by crystallite_subdt
! based on crystallite_subF0,.._subFp0,.._subLp0
! constitutive_state is internally interpolated with .._subState0
! to account for substepping within _integrateStress
! results in crystallite_Fp,.._Lp
!$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)) & ! all undone crystallites
crystallite_onTrack(g,i,e) = crystallite_integrateStress(g,i,e)
enddo
enddo
enddo
!$OMP END PARALLEL DO
write(6,'(i2,x,a10,x,16(l,x))') crystallite_Niteration,'cryst_onT',crystallite_onTrack
write (6,'(a,/,3(3(f12.7,x)/))') 'Lp of 1 8 1',crystallite_Lp(1:3,:,1,8,1)
! --+>> state integration <<+--
!
! incrementing by crystallite_subdt
! based on constitutive_subState0
! results in constitutive_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)) then ! all undone crystallites
crystallite_converged(g,i,e) = crystallite_updateState(g,i,e)
if (crystallite_converged(g,i,e)) then
!$OMP CRITICAL (distributionState)
debug_StateLoopDistribution(crystallite_Niteration) = &
debug_StateLoopDistribution(crystallite_Niteration) + 1
!$OMP END CRITICAL (distributionState)
endif
endif
enddo
enddo
enddo
!$OMP END PARALLEL DO
enddo ! crystallite convergence loop
write(6,*)
write(6,'(a10,x,16(f6.4,x))') 'cryst_frac',crystallite_subFrac
write(6,'(a10,x,16(f6.4,x))') 'cryst_step',crystallite_subStep
write(6,'(a10,x,16(l,x))') 'cryst_req',crystallite_requested
write(6,'(a10,x,16(l,x))') 'cryst_onT',crystallite_onTrack
write(6,'(a10,x,16(l,x))') 'cryst_cvg',crystallite_converged
write(6,'(a10,x,16(e8.3,x))') 'cryst_dt',crystallite_subdt
enddo ! cutback loop
! ------ check for non-converged crystallites ------
!$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 (.not. crystallite_converged(g,i,e)) then ! respond fully elastically
call IO_warning(600,e,i,g)
invFp = math_inv3x3(crystallite_partionedFp0(:,:,g,i,e))
Fe_guess = math_mul33x33(crystallite_partionedF(:,:,g,i,e),invFp)
PK2 = math_Mandel6to33( &
math_mul66x6( &
0.5_pReal*constitutive_homogenizedC(g,i,e), &
math_Mandel33to6( &
math_mul33x33(transpose(Fe_guess),Fe_guess) - math_I3 &
) &
) &
)
crystallite_P(:,:,g,i,e) = math_mul33x33(Fe_guess,math_mul33x33(PK2,transpose(invFp)))
endif
enddo
enddo
enddo
!$OMP END PARALLEL DO
! ------ stiffness calculation ------
if(updateJaco) then ! Jacobian required
if (debugger) then
!$OMP CRITICAL (write2out)
write (6,*) 'Jacobian calc'
write(6,'(a10,x,16(f6.4,x))') 'cryst_dt',crystallite_subdt
!$OMP END CRITICAL (write2out)
endif
!$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_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...
myF = crystallite_subF(:,:,g,i,e) ! ... and kinematics
myFp = crystallite_Fp(:,:,g,i,e)
myFe = crystallite_Fe(:,:,g,i,e)
myLp = crystallite_Lp(:,:,g,i,e)
myP = crystallite_P(:,:,g,i,e)
do k = 1,3 ! perturbation...
do l = 1,3 ! ...components
crystallite_subF(:,:,g,i,e) = myF ! initialize perturbed F to match converged
crystallite_subF(k,l,g,i,e) = crystallite_subF(k,l,g,i,e) + pert_Fg ! perturb single component
doneAndHappy = .false.
crystallite_Niteration = 0_pInt
do while(.not. doneAndHappy(1) .and. crystallite_Niteration < nCryst) ! keep cycling until done (though unhappy)
crystallite_Niteration = crystallite_Niteration + 1_pInt
doneAndHappy = crystallite_integrateStress(g,i,e) ! stress of perturbed situation (overwrites _P,_Tstar_v,_Fp,_Lp,_Fe)
if (doneAndHappy(2)) & ! happy stress allows for state update
doneAndHappy = crystallite_updateState(g,i,e)
end do
if (doneAndHappy(2)) & ! happy outcome warrants stiffness update
crystallite_dPdF(:,:,k,l,g,i,e) = (crystallite_p(:,:,g,i,e) - myP)/pert_Fg ! tangent dP_ij/dFg_kl
!$OMP CRITICAL (out)
debug_StiffnessStateLoopDistribution(crystallite_Niteration) = &
debug_StiffnessstateLoopDistribution(crystallite_Niteration) + 1
!$OMP END CRITICAL (out)
end do
end do
constitutive_state(g,i,e)%p = myState ! restore unperturbed, converged state...
crystallite_Fp(:,:,g,i,e) = myFp ! ... and kinematics
crystallite_Fe(:,:,g,i,e) = myFe
crystallite_Lp(:,:,g,i,e) = myLp
crystallite_P(:,:,g,i,e) = myP
else ! grain has not converged
crystallite_dPdF(:,:,:,:,g,i,e) = crystallite_fallbackdPdF(:,:,:,:,g,i,e) ! use fallback
endif
enddo
enddo
enddo
!$OMP END PARALLEL DO
endif
end subroutine
!********************************************************************
! update the internal state of the constitutive law
! and tell whether state has converged
!********************************************************************
function crystallite_updateState(&
g,& ! grain number
i,& ! integration point number
e & ! element number
)
use prec, only: pReal,pInt,rTol_crystalliteState
use constitutive, only: constitutive_dotState,constitutive_sizeDotState,&
constitutive_subState0,constitutive_state
use debug
logical crystallite_updateState
integer(pLongInt) tick,tock,tickrate,maxticks
integer(pInt) g,i,e,mySize
real(pReal), dimension(6) :: Tstar_v
real(pReal), dimension(constitutive_sizeDotState(g,i,e)) :: residuum
mySize = constitutive_sizeDotState(g,i,e)
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) ! residuum from evolution of microstructure
call system_clock(count=tock,count_rate=tickrate,count_max=maxticks)
debug_cumDotStateCalls = debug_cumDotStateCalls + 1_pInt
debug_cumDotStateTicks = debug_cumDotStateTicks + tock-tick
if (tock < tick) debug_cumDotStateTicks = debug_cumDotStateTicks + maxticks
constitutive_state(g,i,e)%p(1:mySize) = constitutive_state(g,i,e)%p(1:mySize) - residuum ! update of microstructure
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
return
end function
!***********************************************************************
!*** calculation of stress (P), stiffness (dPdF), ***
!*** and announcement of any ***
!*** acceleration of the Newton-Raphson correction ***
!***********************************************************************
function crystallite_integrateStress(&
g,& ! grain number
i,& ! integration point number
e) ! element number
use prec, only: pReal,pInt,subStepMin
use debug
use constitutive, only: constitutive_state,constitutive_subState0,constitutive_sizeState
use math
implicit none
integer(pInt), intent(in) :: e,i,g
integer(pInt) mySize, Niteration
real(pReal) dt_aim,subFrac,subStep,det
real(pReal), dimension(3,3) :: inv
real(pReal), dimension(3,3) :: Fg_current,Fg_aim,deltaFg
real(pReal), dimension(3,3) :: Fp_current,Fp_new
real(pReal), dimension(3,3) :: Fe_current,Fe_new
real(pReal), dimension(constitutive_sizeState(g,i,e)) :: interpolatedState
logical crystallite_integrateStress ! still on track?
mySize = constitutive_sizeState(g,i,e)
deltaFg = crystallite_subF(:,:,g,i,e) - crystallite_subF0(:,:,g,i,e)
Fg_current = crystallite_subF0(:,:,g,i,e) ! initialize to start of inc
Fp_current = crystallite_subFp0(:,:,g,i,e)
Fe_current = math_mul33x33(Fg_current,math_inv3x3(Fp_current))
subFrac = 0.0_pReal
subStep = 1.0_pReal
Niteration = 0_pInt
crystallite_integrateStress = .false. ! be pessimisitc
! begin the cutback loop
do while (subStep > subStepMin) ! continue until finished or too much cut backing
Niteration = Niteration + 1
Fg_aim = Fg_current + subStep*deltaFg ! aim for Fg
dt_aim = subStep*crystallite_subdt(g,i,e) ! aim for dt
debugger = (g == 1 .and. i == 8 .and. e == 1)
call TimeIntegration(crystallite_integrateStress,&
crystallite_Lp(:,:,g,i,e),crystallite_Fp(:,:,g,i,e),crystallite_Fe(:,:,g,i,e),&
crystallite_Tstar_v(:,g,i,e),crystallite_P(:,:,g,i,e), &
Fg_aim,Fp_current,crystallite_Temperature(g,i,e),&
( subFrac+subStep)*constitutive_state(g,i,e)%p(1:mySize) + &
(1.0_pReal-subFrac-subStep)*constitutive_subState0(g,i,e)%p(1:mySize),& ! interpolated state
dt_aim,g,i,e)
if (crystallite_integrateStress) then ! happy with time integration
if (e == 1 .and. i == 8 .and. g == 1) then
write(6,*) '*** winding forward in IntegrateStress ***'
write(6,*) subFrac, subFrac+subStep
write (6,'(a,/,3(3(f12.7,x)/))') 'Lp of 1 8 1',crystallite_Lp(1:3,:,1,8,1)
write (6,'(a,/,3(3(f12.7,x)/))') 'Fp of 1 8 1',crystallite_Fp(1:3,:,1,8,1)
endif
Fg_current = Fg_aim ! wind forward
Fe_current = crystallite_Fe(:,:,g,i,e)
Fp_current = crystallite_Fp(:,:,g,i,e)
subFrac = subFrac + subStep
subStep = min(1.0_pReal-subFrac, subStep*2.0_pReal) ! accelerate
else ! time integration encountered trouble
subStep = 0.5_pReal * subStep ! cut time step in half
crystallite_Lp(:,:,g,i,e) = 0.5_pReal*(crystallite_Lp(:,:,g,i,e) + & ! interpolate Lp and L
(math_I3 - math_mul33x33(Fp_current,&
math_mul33x33(math_inv3x3(Fg_aim),Fe_current)))/dt_aim)
endif
enddo ! potential substepping
!$OMP CRITICAL (distributionStress)
debug_StressLoopDistribution(Niteration) = debug_StressLoopDistribution(Niteration) + 1
!$OMP END CRITICAL (distributionStress)
return ! with final happyness
end function
!***********************************************************************
!*** fully-implicit two-level time integration ***
!*** based on a residuum in Lp and intermediate ***
!*** acceleration of the Newton-Raphson correction ***
!***********************************************************************
subroutine TimeIntegration(&
happy,& ! return status
Lpguess,& ! guess of plastic velocity gradient
Fp_new,& ! new plastic deformation gradient
Fe_new,& ! new "elastic" deformation gradient
Tstar_v,& ! Stress vector
P,& ! 1st PK stress (taken as initial guess if /= 0)
!
Fg_new,& ! new total def gradient
Fp_old,& ! former plastic def gradient
Temperature,& ! temperature
state,& ! microstructural state
dt,& ! time increment
grain,& ! grain number
ip,& ! integration point number
cp_en & ! element number
)
use prec
use debug
use mesh, only: mesh_element
use constitutive, only: constitutive_microstructure,constitutive_homogenizedC,constitutive_LpAndItsTangent,&
constitutive_sizeState
use math
use IO
implicit none
logical, intent(out) :: happy
real(pReal), dimension(3,3), intent(inout) :: Lpguess
real(pReal), dimension(3,3), intent(out) :: Fp_new,Fe_new,P
real(pReal), dimension(6), intent(out) :: Tstar_v
real(pReal), dimension(3,3), intent(in) :: Fg_new,Fp_old
real(pReal), intent(in) :: Temperature,dt
integer(pInt), intent(in) :: cp_en, ip, grain
real(pReal), dimension(constitutive_sizeState(grain,ip,cp_en)), intent(in) :: state
logical error
integer(pInt) Niteration,dummy, i,j,k,l,m,n
integer(pLongInt) tick,tock,tickrate,maxticks
real(pReal) p_hydro,det, leapfrog,maxleap
real(pReal), dimension(3,3,3,3) :: C
real(pReal), dimension(9,9) :: dLp,dTdLp,dRdLp,invdRdLp,eye2
real(pReal), dimension(6,6) :: C_66
real(pReal), dimension(3,3) :: invFp_new,invFp_old,Lp,Lpguess_old,Rinner,Rinner_old,A,B,BT,AB,BTA
happy = .false. ! be pessimistic
eye2 = math_identity2nd(9)
invFp_old = math_inv3x3(Fp_old) ! inversion of Fp_old
if (all(invFp_old == 0.0_pReal)) return ! failed
! write (6,'(a,/,3(3(f12.7,x)/))') 'Fp old',Fp_old
! write (6,'(a,/,3(3(f12.7,x)/))') 'Fp old inv',invFp_old
A = math_mul33x33(transpose(invFp_old), math_mul33x33(transpose(Fg_new),math_mul33x33(Fg_new,invFp_old)))
!$OMP CRITICAL (write2out)
! if (debugger) write (6,'(a,/,3(3(f12.7,x)/))') 'Fg to be calculated',Fg_new
!$OMP END CRITICAL (write2out)
call constitutive_microstructure(Temperature,grain,ip,cp_en)
C_66 = constitutive_homogenizedC(grain,ip,cp_en)
C = math_Mandel66to3333(C_66) ! 4th rank elasticity tensor
Niteration = 0_pInt
leapfrog = 1.0_pReal ! correction as suggested by invdRdLp-step
maxleap = 1024.0_pReal ! preassign maximum acceleration level
Lpguess_old = Lpguess ! consider present Lpguess good (i.e. worth remembering)
Inner: do ! inner iteration: Lp
Niteration = Niteration + 1
if (Niteration > nLp) then ! too many loops required
Lpguess = Lpguess_old ! do not trust the last update but resort to former one
return
endif
! write(6,*) 'iteration',Niteration
! write (6,'(a,/,3(3(f12.7,x)/))') 'Lpguess',Lpguess
B = math_i3 - dt*Lpguess
BT = transpose(B)
AB = math_mul33x33(A,B)
BTA = math_mul33x33(BT,A)
Tstar_v = 0.5_pReal*math_mul66x6(C_66,math_mandel33to6(math_mul33x33(BT,AB)-math_I3))
p_hydro = (Tstar_v(1) + Tstar_v(2) + Tstar_v(3))/3.0_pReal
forall(i=1:3) Tstar_v(i) = Tstar_v(i) - p_hydro ! subtract hydrostatic pressure
! write (6,'(a,/,6(f12.7,x))') 'Tstar',Tstar_v
call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
call constitutive_LpAndItsTangent(Lp,dLp, Tstar_v,Temperature,grain,ip,cp_en)
call system_clock(count=tock,count_rate=tickrate,count_max=maxticks)
debug_cumLpCalls = debug_cumLpCalls + 1_pInt
debug_cumLpTicks = debug_cumLpTicks + tock-tick
if (tock < tick) debug_cumLpTicks = debug_cumLpTicks + maxticks
Rinner = Lpguess - Lp ! update current residuum
! write (6,'(a,/,3(3(f12.7,x)/))') 'Lp',Lp
! write (6,'(a,/,3(3(f12.7,x)/))') 'Residuum',Rinner
if (.not.(any(Rinner/=Rinner)) .and. & ! exclude any NaN in residuum
( ( maxval(abs(Rinner)) < aTol_crystalliteStress) .or. & ! below abs tol .or.
( any(abs(dt*Lpguess) > relevantStrain) .and. & ! worth checking? .and.
maxval(abs(Rinner/Lpguess),abs(dt*Lpguess) > relevantStrain) < rTol_crystalliteStress & ! below rel tol
) &
) &
) &
exit Inner ! convergence
!
! check for acceleration/deceleration in Newton--Raphson correction
!
if (any(Rinner/=Rinner) .and. & ! NaN occured at regular speed
leapfrog == 1.0) then
Lpguess = Lpguess_old ! restore known good guess and croak for cutback
return
elseif (leapfrog > 1.0_pReal .and. & ! at fast pace ?
(sum(Rinner*Rinner) > sum(Rinner_old*Rinner_old) .or. & ! worse residuum
sum(Rinner*Rinner_old) < 0.0_pReal) .or. & ! residuum changed sign (overshoot)
any(Rinner/=Rinner) ) then ! NaN
maxleap = 0.5_pReal * leapfrog ! limit next acceleration
leapfrog = 1.0_pReal ! grinding halt
else ! better residuum
dTdLp = 0.0_pReal ! calc dT/dLp
forall (i=1:3,j=1:3,k=1:3,l=1:3,m=1:3,n=1:3) &
dTdLp(3*(i-1)+j,3*(k-1)+l) = dTdLp(3*(i-1)+j,3*(k-1)+l) + &
C(i,j,l,n)*AB(k,n)+C(i,j,m,l)*BTA(m,k)
dTdLp = -0.5_pReal*dt*dTdLp
dRdLp = eye2 - math_mul99x99(dLp,dTdLp) ! calc dR/dLp
invdRdLp = 0.0_pReal
call math_invert(9,dRdLp,invdRdLp,dummy,error) ! invert dR/dLp --> dLp/dR
if (error) then
if (debugger) then
!$OMP CRITICAL (write2out)
write (6,*) 'inversion dR/dLp failed',grain,ip,cp_en
! write (6,'(a,/,9(9(e9.3,x)/))') 'dRdLp', dRdLp(1:9,:)
! write (6,'(a,/,3(4(f9.3,x)/))') 'state_new / MPa',state/1e6_pReal
write (6,'(a,/,3(3(f12.7,x)/))') 'Lpguess',Lpguess(1:3,:)
write (6,'(a,/,3(3(e12.7,x)/))') 'Lp',Lp(1:3,:)
write (6,'(a,/,6(f9.3,x))') 'Tstar / MPa',Tstar_v/1e6_pReal
!$OMP END CRITICAL (write2out)
endif
return
endif
Rinner_old = Rinner ! remember current residuum
Lpguess_old = Lpguess ! remember current Lp guess
if (Niteration > 1 .and. leapfrog < maxleap) leapfrog = 2.0_pReal * leapfrog ! accelerate if ok
endif
Lpguess = Lpguess_old ! start from current guess
Rinner = Rinner_old ! use current residuum
forall (i=1:3,j=1:3,k=1:3,l=1:3) & ! leapfrog to updated Lpguess
Lpguess(i,j) = Lpguess(i,j) - leapfrog*invdRdLp(3*(i-1)+j,3*(k-1)+l)*Rinner(k,l)
enddo Inner
!$OMP CRITICAL (distributionLp)
debug_LpLoopDistribution(Niteration) = debug_LpLoopDistribution(Niteration) + 1
!$OMP END CRITICAL (distributionLp)
invFp_new = math_mul33x33(invFp_old,B)
if (debugger) then
write (6,'(a,x,f10.6,/,3(3(f12.7,x)/))') 'Lp(guess)',dt,Lpguess(1:3,:)
write (6,'(a,/,3(3(f12.7,x)/))') 'invFp_old',invFp_old(1:3,:)
write (6,'(a,/,3(3(f12.7,x)/))') 'B',B(1:3,:)
write (6,'(a,/,3(3(f12.7,x)/))') 'invFp_new',invFp_new(1:3,:)
endif
call math_invert3x3(invFp_new,Fp_new,det,error)
if (debugger) then
write (6,'(a,/,3(3(f12.7,x)/))') 'invFp_new',invFp_new(1:3,:)
write (6,'(a,/,3(3(f12.7,x)/))') 'Fp_new',Fp_new(1:3,:)
write (6,'(a,x,l,x,a,f10.6)') 'with inversion error:',error,'and determinant:',det
endif
if (error) return ! inversion failed
Fp_new = Fp_new*det**(1.0_pReal/3.0_pReal) ! regularize Fp by det = det(InvFp_new) !!
forall (i=1:3) Tstar_v(i) = Tstar_v(i) + p_hydro ! add hydrostatic component back
Fe_new = math_mul33x33(Fg_new,invFp_new) ! calc resulting Fe
P = math_mul33x33(Fe_new,math_mul33x33(math_Mandel6to33(Tstar_v),transpose(invFp_new))) ! first PK stress
happy = .true. ! now smile...
return
end subroutine
!********************************************************************
! return results of particular grain
!********************************************************************
function crystallite_postResults(&
Tstar_v,& ! stress
Temperature, & ! temperature
dt,& ! time increment
g,& ! grain number
i,& ! integration point number
e & ! element number
)
use prec, only: pInt,pReal
use math, only: math_pDecomposition,math_RtoEuler, inDeg
use IO, only: IO_warning
use material, only: material_phase,material_volfrac
use constitutive, only: constitutive_sizePostResults, constitutive_postResults
implicit none
integer(pInt), intent(in) :: g,i,e
real(pReal), intent(in) :: Temperature,dt
real(pReal), dimension(6), intent(in) :: Tstar_v
real(pReal), dimension(3,3) :: U,R
logical error
real(pReal), dimension(crystallite_Nresults + constitutive_sizePostResults(g,i,e)) :: crystallite_postResults
if (crystallite_Nresults >= 2) then
crystallite_postResults(1) = material_phase(g,i,e)
crystallite_postResults(2) = material_volfrac(g,i,e)
endif
if (crystallite_Nresults >= 5) then
call math_pDecomposition(crystallite_Fe(:,:,g,i,e),U,R,error) ! polar decomposition of Fe
if (error) then
call IO_warning(650,e,i,g)
crystallite_postResults(3:5) = (/400.0,400.0,400.0/) ! fake orientation
else
crystallite_postResults(3:5) = math_RtoEuler(transpose(R))*inDeg ! orientation
endif
endif
crystallite_postResults(crystallite_Nresults+1:crystallite_Nresults+constitutive_sizePostResults(g,i,e)) = &
constitutive_postResults(Tstar_v,Temperature,dt,g,i,e)
return
end function
END MODULE
!##############################################################

81
trunk/debug.f90
View File

@ -6,9 +6,10 @@
implicit none implicit none
integer(pInt), dimension(nCutback+1) :: debug_cutbackDistribution = 0_pInt integer(pInt), dimension(nLp) :: debug_LpLoopDistribution = 0_pInt
integer(pInt), dimension(nInner) :: debug_InnerLoopDistribution = 0_pInt integer(pInt), dimension(nStress) :: debug_StressLoopDistribution = 0_pInt
integer(pInt), dimension(nOuter) :: debug_OuterLoopDistribution = 0_pInt integer(pInt), dimension(nCryst) :: debug_StateLoopDistribution = 0_pInt
integer(pInt), dimension(nCryst) :: debug_StiffnessStateLoopDistribution = 0_pInt
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(pInt) :: debug_cumLpCalls = 0_pInt integer(pInt) :: debug_cumLpCalls = 0_pInt
@ -18,6 +19,24 @@
CONTAINS CONTAINS
!********************************************************************
! reset debug distributions
!********************************************************************
SUBROUTINE debug_reset()
use prec
implicit none
debug_LpLoopDistribution = 0_pInt ! initialize debugging data
debug_StressLoopDistribution = 0_pInt
debug_StateLoopDistribution = 0_pInt
debug_StiffnessStateLoopDistribution = 0_pInt
debug_cumLpTicks = 0_pInt
debug_cumDotStateTicks = 0_pInt
debug_cumLpCalls = 0_pInt
debug_cumDotStateCalls = 0_pInt
END SUBROUTINE
!******************************************************************** !********************************************************************
! write debug statements to standard out ! write debug statements to standard out
@ -47,34 +66,50 @@
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,*) 'distribution_cutback :'
do i=0,nCutback
if (debug_cutbackDistribution(i+1) /= 0) write(6,*) i,debug_cutbackDistribution(i+1)
enddo
write(6,*) 'total',sum(debug_cutbackDistribution)
write(6,*)
integral = 0_pInt integral = 0_pInt
write(6,*) 'distribution_InnerLoop :' write(6,*)
do i=1,nInner write(6,*) 'distribution_LpLoop :'
if (debug_InnerLoopDistribution(i) /= 0) then do i=1,nLp
integral = integral + i*debug_InnerLoopDistribution(i) if (debug_LpLoopDistribution(i) /= 0) then
write(6,*) i,debug_InnerLoopDistribution(i) integral = integral + i*debug_LpLoopDistribution(i)
write(6,*) i,debug_LpLoopDistribution(i)
endif endif
enddo enddo
write(6,*) 'total',sum(debug_InnerLoopDistribution),integral write(6,*) 'total',sum(debug_LpLoopDistribution),integral
write(6,*)
integral = 0_pInt integral = 0_pInt
write(6,*) 'distribution_OuterLoop :' write(6,*)
do i=1,nOuter write(6,*) 'distribution_StressLoop :'
if (debug_OuterLoopDistribution(i) /= 0) then do i=1,nStress
integral = integral + i*debug_OuterLoopDistribution(i) if (debug_StressLoopDistribution(i) /= 0) then
write(6,*) i,debug_OuterLoopDistribution(i) integral = integral + i*debug_StressLoopDistribution(i)
write(6,*) i,debug_StressLoopDistribution(i)
endif endif
enddo enddo
write(6,*) 'total',sum(debug_OuterLoopDistribution),integral write(6,*) 'total',sum(debug_StressLoopDistribution),integral
integral = 0_pInt
write(6,*)
write(6,*) 'distribution_StateLoop :'
do i=1,nCryst
if (debug_StateLoopDistribution(i) /= 0) then
integral = integral + i*debug_StateLoopDistribution(i)
write(6,*) i,debug_StateLoopDistribution(i)
endif
enddo
write(6,*) 'total',sum(debug_StateLoopDistribution),integral
integral = 0_pInt
write(6,*)
write(6,*) 'distribution_StiffnessStateLoop :'
do i=1,nCryst
if (debug_StiffnessStateLoopDistribution(i) /= 0) then
integral = integral + i*debug_StiffnessStateLoopDistribution(i)
write(6,*) i,debug_StiffnessStateLoopDistribution(i)
endif
enddo
write(6,*) 'total',sum(debug_StiffnessStateLoopDistribution),integral
write(6,*) write(6,*)
END SUBROUTINE END SUBROUTINE

505
trunk/homogenization.f90 Normal file
View File

@ -0,0 +1,505 @@
!***************************************
!* Module: HOMOGENIZATION *
!***************************************
!* contains: *
!* - _init *
!* - materialpoint_stressAndItsTangent *
!* - _partitionDeformation *
!* - _updateState *
!* - _averageStressAndItsTangent *
!* - _postResults *
!***************************************
MODULE homogenization
!*** Include other modules ***
use prec, only: pInt,pReal,p_vec
implicit none
! ****************************************************************
! *** General variables for the homogenization at a ***
! *** material point ***
! ****************************************************************
type(p_vec), dimension(:,:), allocatable :: homogenization_state0, & ! pointer array to homogenization state at start of FE increment
homogenization_subState0, & ! pointer array to homogenization state at start of homogenization increment
homogenization_state ! pointer array to current homogenization state (end of converged time step)
integer(pInt), dimension(:,:), allocatable :: homogenization_sizeState, & ! size of state array per grain
homogenization_sizePostResults ! size of postResults array per material point
real(pReal), dimension(:,:,:,:,:,:), allocatable :: materialpoint_dPdF ! tangent of first P--K stress at IP
real(pReal), dimension(:,:,:,:), allocatable :: materialpoint_F0, & ! def grad of IP at start of FE increment
materialpoint_F, & ! def grad of IP to be reached at end of FE increment
materialpoint_subF0, & ! def grad of IP at beginning of homogenization increment
materialpoint_subF, & ! def grad of IP to be reached at end of homog inc
materialpoint_P ! first P--K stress of IP
real(pReal), dimension(:,:), allocatable :: materialpoint_Temperature, & ! temperature at IP
materialpoint_subFrac, &
materialpoint_subStep, &
materialpoint_subdt
real(pReal), dimension(:,:,:), allocatable :: materialpoint_results ! results array of material point
logical, dimension(:,:), allocatable :: materialpoint_requested, &
materialpoint_converged
logical, dimension(:,:,:), allocatable :: materialpoint_doneAndHappy
integer(pInt) homogenization_maxSizeState,homogenization_maxSizePostResults
CONTAINS
!**************************************
!* Module initialization *
!**************************************
subroutine homogenization_init()
use prec, only: pReal,pInt
use math, only: math_I3
use IO, only: IO_error, IO_open_file
use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips
use material
use constitutive, only: constitutive_maxSizePostResults
use crystallite, only: crystallite_Nresults
use homogenization_isostrain
! use homogenization_RGC
integer(pInt), parameter :: fileunit = 200
integer(pInt) e,i,g,myInstance
if(.not. IO_open_file(fileunit,material_configFile)) call IO_error (100) ! corrupt config file
call homogenization_isostrain_init(fileunit) ! parse all homogenizations of this type
close(fileunit)
allocate(homogenization_state0(mesh_maxNips,mesh_NcpElems))
allocate(homogenization_subState0(mesh_maxNips,mesh_NcpElems))
allocate(homogenization_state(mesh_maxNips,mesh_NcpElems))
allocate(homogenization_sizeState(mesh_maxNips,mesh_NcpElems)); homogenization_sizeState = 0_pInt
allocate(homogenization_sizePostResults(mesh_maxNips,mesh_NcpElems)); homogenization_sizePostResults = 0_pInt
allocate(materialpoint_dPdF(3,3,3,3,mesh_maxNips,mesh_NcpElems)); materialpoint_dPdF = 0.0_pReal
allocate(materialpoint_F0(3,3,mesh_maxNips,mesh_NcpElems));
allocate(materialpoint_F(3,3,mesh_maxNips,mesh_NcpElems)); materialpoint_F = 0.0_pReal
allocate(materialpoint_subF0(3,3,mesh_maxNips,mesh_NcpElems)); materialpoint_subF0 = 0.0_pReal
allocate(materialpoint_subF(3,3,mesh_maxNips,mesh_NcpElems)); materialpoint_subF = 0.0_pReal
allocate(materialpoint_P(3,3,mesh_maxNips,mesh_NcpElems)); materialpoint_P = 0.0_pReal
allocate(materialpoint_Temperature(mesh_maxNips,mesh_NcpElems)); materialpoint_Temperature = 0.0_pReal
allocate(materialpoint_subFrac(mesh_maxNips,mesh_NcpElems)); materialpoint_subFrac = 0.0_pReal
allocate(materialpoint_subStep(mesh_maxNips,mesh_NcpElems)); materialpoint_subStep = 0.0_pReal
allocate(materialpoint_subdt(mesh_maxNips,mesh_NcpElems)); materialpoint_subdt = 0.0_pReal
allocate(materialpoint_requested(mesh_maxNips,mesh_NcpElems)); materialpoint_requested = .false.
allocate(materialpoint_converged(mesh_maxNips,mesh_NcpElems)); materialpoint_converged = .true.
allocate(materialpoint_doneAndHappy(2,mesh_maxNips,mesh_NcpElems)); materialpoint_doneAndHappy = .true.
forall (i = 1:mesh_maxNips,e = 1:mesh_NcpElems)
materialpoint_F0(:,:,i,e) = math_I3
materialpoint_F(:,:,i,e) = math_I3
end forall
do e = 1,mesh_NcpElems ! loop over elements
myInstance = homogenization_typeInstance(mesh_element(3,e))
do i = 1,FE_Nips(mesh_element(2,e)) ! loop over IPs
select case(homogenization_type(mesh_element(3,e)))
case (homogenization_isostrain_label)
if (homogenization_isostrain_sizeState(myInstance) > 0_pInt) then
allocate(homogenization_state0(i,e)%p(homogenization_isostrain_sizeState(myInstance)))
allocate(homogenization_subState0(i,e)%p(homogenization_isostrain_sizeState(myInstance)))
allocate(homogenization_state(i,e)%p(homogenization_isostrain_sizeState(myInstance)))
homogenization_state0(i,e)%p = homogenization_isostrain_stateInit(myInstance)
homogenization_sizeState(i,e) = homogenization_isostrain_sizeState(myInstance)
endif
homogenization_sizePostResults(i,e) = homogenization_isostrain_sizePostResults(myInstance)
case default
call IO_error(200,ext_msg=homogenization_type(mesh_element(3,e))) ! unknown type 200 is phase!
end select
enddo
enddo
homogenization_maxSizeState = maxval(homogenization_sizeState)
homogenization_maxSizePostResults = maxval(homogenization_sizePostResults)
allocate(materialpoint_results( 1+homogenization_maxSizePostResults + &
homogenization_maxNgrains*(1+crystallite_Nresults+constitutive_maxSizePostResults), mesh_maxNips,mesh_NcpElems))
! *** Output to MARC output file ***
!$OMP CRITICAL (write2out)
write(6,*)
write(6,*) '<<<+- homogenization init -+>>>'
write(6,*)
write(6,'(a32,x,7(i5,x))') 'homogenization_state0: ', shape(homogenization_state0)
write(6,'(a32,x,7(i5,x))') 'homogenization_subState0: ', shape(homogenization_subState0)
write(6,'(a32,x,7(i5,x))') 'homogenization_state: ', shape(homogenization_state)
write(6,'(a32,x,7(i5,x))') 'homogenization_sizeState: ', shape(homogenization_sizeState)
write(6,'(a32,x,7(i5,x))') 'homogenization_sizePostResults: ', shape(homogenization_sizePostResults)
write(6,*)
write(6,'(a32,x,7(i5,x))') 'materialpoint_dPdF: ', shape(materialpoint_dPdF)
write(6,'(a32,x,7(i5,x))') 'materialpoint_F0: ', shape(materialpoint_F0)
write(6,'(a32,x,7(i5,x))') 'materialpoint_F: ', shape(materialpoint_F)
write(6,'(a32,x,7(i5,x))') 'materialpoint_subF0: ', shape(materialpoint_subF0)
write(6,'(a32,x,7(i5,x))') 'materialpoint_subF: ', shape(materialpoint_subF)
write(6,'(a32,x,7(i5,x))') 'materialpoint_P: ', shape(materialpoint_P)
write(6,'(a32,x,7(i5,x))') 'materialpoint_Temperature: ', shape(materialpoint_Temperature)
write(6,'(a32,x,7(i5,x))') 'materialpoint_subFrac: ', shape(materialpoint_subFrac)
write(6,'(a32,x,7(i5,x))') 'materialpoint_subStep: ', shape(materialpoint_subStep)
write(6,'(a32,x,7(i5,x))') 'materialpoint_subdt: ', shape(materialpoint_subdt)
write(6,'(a32,x,7(i5,x))') 'materialpoint_requested: ', shape(materialpoint_requested)
write(6,'(a32,x,7(i5,x))') 'materialpoint_converged: ', shape(materialpoint_converged)
write(6,'(a32,x,7(i5,x))') 'materialpoint_doneAndHappy: ', shape(materialpoint_doneAndHappy)
write(6,*)
write(6,'(a32,x,7(i5,x))') 'materialpoint_results: ', shape(materialpoint_results)
write(6,*)
write(6,'(a32,x,7(i5,x))') 'maxSizeState: ', homogenization_maxSizeState
write(6,'(a32,x,7(i5,x))') 'maxSizePostResults: ', homogenization_maxSizePostResults
call flush(6)
!$OMP END CRITICAL (write2out)
return
end subroutine
!********************************************************************
!* parallelized calculation of
!* stress and corresponding tangent
!* at material points
!********************************************************************
subroutine materialpoint_stressAndItsTangent(&
updateJaco,& ! flag to initiate Jacobian updating
dt & ! time increment
)
use prec, only: pInt,pReal, subStepMin,nHomog
use FEsolving, only: FEsolving_execElem, FEsolving_execIP
use mesh, only: mesh_element
use material, only: homogenization_Ngrains
use constitutive, only: constitutive_state0, constitutive_partionedState0, constitutive_state
use crystallite
implicit none
real(pReal), intent(in) :: dt
logical, intent(in) :: updateJaco
integer(pInt) homogenization_Niteration
integer(pInt) g,i,e,myNgrains
! ------ initialize to starting condition ------
write (6,*)
write (6,*) 'Material Point start'
write (6,'(a,/,3(3(f12.7,x)/))') 'F0 of 8 1',materialpoint_F0(1:3,:,8,1)
write (6,'(a,/,3(3(f12.7,x)/))') 'F of 8 1',materialpoint_F(1:3,:,8,1)
write (6,'(a,/,3(3(f12.7,x)/))') 'Fp0 of 1 8 1',crystallite_Fp0(1:3,:,1,8,1)
write (6,'(a,/,3(3(f12.7,x)/))') 'Lp0 of 1 8 1',crystallite_Lp0(1:3,:,1,8,1)
!$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
! initialize restoration points of grain...
forall (g = 1:myNgrains) constitutive_partionedState0(g,i,e)%p = constitutive_state0(g,i,e)%p ! ...microstructures
crystallite_partionedFp0(:,:,1:myNgrains,i,e) = crystallite_Fp0(:,:,1:myNgrains,i,e) ! ...plastic def grads
crystallite_partionedLp0(:,:,1:myNgrains,i,e) = crystallite_Lp0(:,:,1:myNgrains,i,e) ! ...plastic velocity grads
crystallite_partionedF0(:,:,1:myNgrains,i,e) = crystallite_F0(:,:,1:myNgrains,i,e) ! ...def grads
! initialize restoration points of ...
if (homogenization_sizeState(i,e) > 0_pInt) &
homogenization_subState0(i,e)%p = homogenization_state0(i,e)%p ! ...internal homogenizaiton state
materialpoint_subF0(:,:,i,e) = materialpoint_F0(:,:,i,e) ! ...def grad
materialpoint_subFrac(i,e) = 0.0_pReal
materialpoint_subStep(i,e) = 2.0_pReal
materialpoint_converged(i,e) = .false. ! pretend failed step of twice the required size
materialpoint_requested(i,e) = .true. ! everybody requires calculation
enddo
enddo
!$OMP END PARALLEL DO
! ------ cutback loop ------
do while (any(materialpoint_subStep(:,FEsolving_execELem(1):FEsolving_execElem(2)) > subStepMin)) ! cutback loop for material points
!$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
if (materialpoint_converged(i,e)) then
materialpoint_subFrac(i,e) = materialpoint_subFrac(i,e) + materialpoint_subStep(i,e)
materialpoint_subStep(i,e) = min(1.0_pReal-materialpoint_subFrac(i,e), 2.0_pReal * materialpoint_subStep(i,e))
if (materialpoint_subStep(i,e) > subStepMin) then ! still stepping needed
! wind forward grain starting point of...
crystallite_partionedF0(:,:,1:myNgrains,i,e) = crystallite_partionedF(:,:,1:myNgrains,i,e) ! ...def grads
crystallite_partionedFp0(:,:,1:myNgrains,i,e) = crystallite_Fp(:,:,1:myNgrains,i,e) ! ...plastic def grads
crystallite_partionedLp0(:,:,1:myNgrains,i,e) = crystallite_Lp(:,:,1:myNgrains,i,e) ! ...plastic velocity grads
forall (g = 1:myNgrains) constitutive_partionedState0(g,i,e)%p = constitutive_state(g,i,e)%p ! ...microstructures
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
endif
else
materialpoint_subStep(i,e) = 0.5_pReal * materialpoint_subStep(i,e) ! cut step in half and restore...
! ####### why not resetting F0 ?!?!?
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
forall (g = 1:myNgrains) constitutive_state(g,i,e)%p = constitutive_partionedState0(g,i,e)%p ! ...microstructures
if (homogenization_sizeState(i,e) > 0_pInt) &
homogenization_state(i,e)%p = homogenization_subState0(i,e)%p ! ...internal state of homog scheme
endif
materialpoint_requested(i,e) = materialpoint_subStep(i,e) > subStepMin
if (materialpoint_requested(i,e)) then
materialpoint_subF(:,:,i,e) = materialpoint_subF0(:,:,i,e) + &
materialpoint_subStep(i,e) * (materialpoint_F(:,:,i,e) - materialpoint_F0(:,:,i,e))
materialpoint_subdt(i,e) = materialpoint_subStep(i,e) * dt
materialpoint_doneAndHappy(:,i,e) = (/.false.,.true./)
endif
enddo
enddo
!$OMP END PARALLEL DO
! tell what is requested
write(6,'(a14,x,16(l,x))') 'matpnt_req',materialpoint_requested
write(6,'(a14,x,16(l,x))') 'matpnt_don',materialpoint_doneAndHappy(1,:,:)
write(6,'(a14,x,16(l,x))') 'matpnt_hpy',materialpoint_doneAndHappy(1,:,:)
write(6,'(a14,x,16(f6.4,x))') 'matpnt_frac',materialpoint_subFrac
write(6,'(a14,x,16(f6.4,x))') 'matpnt_step',materialpoint_subStep
write(6,'(a10,x,16(e8.3,x))') 'matpnt_dt',materialpoint_subdt
write(6,*)
write (6,'(a,/,3(3(f12.7,x)/))') 'subF0 of 8 1',materialpoint_subF0(1:3,:,8,1)
write (6,'(a,/,3(3(f12.7,x)/))') 'subF of 8 1',materialpoint_subF(1:3,:,8,1)
! ------ convergence loop material point homogenization ------
homogenization_Niteration = 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
! --+>> deformation partitioning <<+--
!
! based on materialpoint_subF0,.._subF,
! crystallite_partionedF0,
! homogenization_state
! results in crystallite_partionedF
!$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
if ( materialpoint_requested(i,e) .and. & ! process requested but...
.not. materialpoint_doneAndHappy(1,i,e)) then ! ...not yet done material points
call homogenization_partitionDeformation(i,e) ! partition deformation onto constituents
crystallite_dt(1:myNgrains,i,e) = materialpoint_subdt(i,e) ! propagate materialpoint dt to grains
crystallite_requested(1:myNgrains,i,e) = .true. ! request calculation for constituents
endif
enddo
enddo
!$OMP END PARALLEL DO
! --+>> crystallite integration <<+--
!
! based on crystallite_partionedF0,.._partionedF
! incrementing by crystallite_dt
call crystallite_stressAndItsTangent(updateJaco) ! request stress and tangent calculation for constituent grains
! --+>> state update <<+--
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
if ( materialpoint_requested(i,e) .and. &
.not. materialpoint_doneAndHappy(1,i,e)) then
materialpoint_doneAndHappy(:,i,e) = homogenization_updateState(i,e)
materialpoint_converged(i,e) = all(materialpoint_doneAndHappy(:,i,e)) ! converged if done and happy
endif
enddo
enddo
!$OMP END PARALLEL DO
enddo ! homogenization convergence loop
enddo ! cutback loop
! check for non-performer: any(.not. converged)
! replace with elastic response ?
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
call homogenization_averageStressAndItsTangent(i,e)
enddo
enddo
!$OMP END PARALLEL DO
write (6,*) 'Material Point finished'
write (6,'(a,/,3(3(f12.7,x)/))') 'Lp of 1 8 1',crystallite_Lp(1:3,:,1,8,1)
! how to deal with stiffness?
return
end subroutine
!********************************************************************
!* parallelized calculation of
!* result array at material points
!********************************************************************
subroutine materialpoint_postResults(dt)
use FEsolving, only: FEsolving_execElem, FEsolving_execIP
use mesh, only: mesh_element
use material, only: homogenization_Ngrains
use constitutive, only: constitutive_sizePostResults, constitutive_postResults
use crystallite
implicit none
real(pReal), intent(in) :: dt
integer(pInt) g,i,e,c,d,myNgrains
!$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
c = 0_pInt
d = homogenization_sizePostResults(i,e)
materialpoint_results(c+1,i,e) = d; c = c+1_pInt ! tell size of homogenization results
if (d > 0_pInt) then ! any homogenization results to mention?
materialpoint_results(c+1:c+d,i,e) = & ! tell homogenization results
homogenization_postResults(i,e); c = c+d
endif
do g = 1,myNgrains !
d = crystallite_Nresults+constitutive_sizePostResults(g,i,e)
materialpoint_results(c+1,i,e) = d; c = c+1_pInt ! tell size of crystallite results
materialpoint_results(c+1:c+d,i,e) = & ! tell crystallite results
crystallite_postResults(crystallite_Tstar_v(:,g,i,e),crystallite_Temperature(g,i,e),dt,g,i,e); c = c+d
enddo
enddo
enddo
!$OMP END PARALLEL DO
end subroutine
!********************************************************************
! partition material point def grad onto constituents
!********************************************************************
subroutine homogenization_partitionDeformation(&
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_partionedF0,crystallite_partionedF
use homogenization_isostrain
implicit none
integer(pInt), intent(in) :: ip,el
select case(homogenization_type(mesh_element(3,el)))
case (homogenization_isostrain_label)
call homogenization_isostrain_partitionDeformation(crystallite_partionedF(:,:,:,ip,el), &
crystallite_partionedF0(:,:,:,ip,el),&
materialpoint_subF(:,:,ip,el),&
homogenization_state(ip,el),ip,el)
end select
end subroutine
!********************************************************************
! update the internal state of the homogenization scheme
! and tell whether "done" and "happy" with result
!********************************************************************
function homogenization_updateState(&
ip, & ! integration point
el & ! element
)
use prec, only: pReal,pInt
use mesh, only: mesh_element
use material, only: homogenization_type, homogenization_maxNgrains
use crystallite, only: crystallite_P,crystallite_dPdF
use homogenization_isostrain
implicit none
integer(pInt), intent(in) :: ip,el
logical, dimension(2) :: homogenization_updateState
select case(homogenization_type(mesh_element(3,el)))
case (homogenization_isostrain_label)
homogenization_updateState = &
homogenization_isostrain_updateState(homogenization_state(ip,el), &
crystallite_P(:,:,:,ip,el),crystallite_dPdF(:,:,:,:,:,ip,el),ip,el)
end select
return
end function
!********************************************************************
! derive average stress and stiffness from constituent quantities
!********************************************************************
subroutine homogenization_averageStressAndItsTangent(&
ip, & ! integration point
el & ! element
)
use prec, only: pReal,pInt
use mesh, only: mesh_element
use material, only: homogenization_type, homogenization_maxNgrains
use crystallite, only: crystallite_P,crystallite_dPdF
use homogenization_isostrain
implicit none
integer(pInt), intent(in) :: ip,el
select case(homogenization_type(mesh_element(3,el)))
case (homogenization_isostrain_label)
call homogenization_isostrain_averageStressAndItsTangent(materialpoint_P(:,:,ip,el), materialpoint_dPdF(:,:,:,:,ip,el),&
crystallite_P(:,:,:,ip,el),crystallite_dPdF(:,:,:,:,:,ip,el),ip,el)
end select
return
end subroutine
!********************************************************************
! return array of homogenization results for post file inclusion
! call only, if homogenization_sizePostResults(ip,el) > 0 !!
!********************************************************************
function homogenization_postResults(&
ip, & ! integration point
el & ! element
)
use prec, only: pReal,pInt
use mesh, only: mesh_element
use material, only: homogenization_type
use homogenization_isostrain
implicit none
!* Definition of variables
integer(pInt), intent(in) :: ip,el
real(pReal), dimension(homogenization_sizePostResults(ip,el)) :: homogenization_postResults
homogenization_postResults = 0.0_pReal
select case (homogenization_type(mesh_element(3,el)))
case (homogenization_isostrain_label)
homogenization_postResults = homogenization_isostrain_postResults(homogenization_state(ip,el),ip,el)
end select
return
end function
END MODULE

275
trunk/homogenization_isostrain.f90 Normal file
View File

@ -0,0 +1,275 @@
!*****************************************************
!* Module: HOMOGENIZATION_ISOSTRAIN *
!*****************************************************
!* contains: *
!*****************************************************
! [isostrain]
! type isostrain
! Ngrains 6
! (output) Ngrains
MODULE homogenization_isostrain
!*** Include other modules ***
use prec, only: pReal,pInt
implicit none
character (len=*), parameter :: homogenization_isostrain_label = 'isostrain'
integer(pInt), dimension(:), allocatable :: homogenization_isostrain_sizeState, &
homogenization_isostrain_sizePostResults, &
homogenization_isostrain_Ngrains
character(len=64), dimension(:,:), allocatable :: homogenization_isostrain_output
CONTAINS
!****************************************
!* - homogenization_isostrain_init
!* - homogenization_isostrain_stateInit
!* - homogenization_isostrain_deformationPartititon
!* - homogenization_isostrain_stateUpdate
!* - homogenization_isostrain_averageStressAndItsTangent
!* - homogenization_isostrain_postResults
!****************************************
!**************************************
!* Module initialization *
!**************************************
subroutine homogenization_isostrain_init(&
file & ! file pointer to material configuration
)
use prec, only: pInt, pReal
use math, only: math_Mandel3333to66, math_Voigt66to3333
use IO
use material
integer(pInt), intent(in) :: file
integer(pInt), parameter :: maxNchunks = 2
integer(pInt), dimension(1+2*maxNchunks) :: positions
integer(pInt) section, maxNinstance, i,j,k,l, output
character(len=64) tag
character(len=1024) line
maxNinstance = count(homogenization_type == homogenization_isostrain_label)
if (maxNinstance == 0) return
allocate(homogenization_isostrain_sizeState(maxNinstance)) ; homogenization_isostrain_sizeState = 0_pInt
allocate(homogenization_isostrain_sizePostResults(maxNinstance)); homogenization_isostrain_sizePostResults = 0_pInt
allocate(homogenization_isostrain_Ngrains(maxNinstance)); homogenization_isostrain_Ngrains = 0_pInt
allocate(homogenization_isostrain_output(maxval(homogenization_Noutput), &
maxNinstance)) ; homogenization_isostrain_output = ''
rewind(file)
line = ''
section = 0
do while (IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization) ! wind forward to <homogenization>
read(file,'(a1024)',END=100) line
enddo
do ! read thru sections of phase part
read(file,'(a1024)',END=100) line
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
if (IO_getTag(line,'[',']') /= '') then ! next section
section = section + 1
output = 0 ! reset output counter
endif
if (section > 0 .and. homogenization_type(section) == homogenization_isostrain_label) then ! one of my sections
i = homogenization_typeInstance(section) ! which instance of my type is present homogenization
positions = IO_stringPos(line,maxNchunks)
tag = IO_lc(IO_stringValue(line,positions,1)) ! extract key
select case(tag)
case ('(output)')
output = output + 1
homogenization_isostrain_output(output,i) = IO_lc(IO_stringValue(line,positions,2))
case ('ngrains')
homogenization_isostrain_Ngrains(i) = IO_intValue(line,positions,2)
end select
endif
enddo
100 do i = 1,maxNinstance ! sanity checks
enddo
do i = 1,maxNinstance
homogenization_isostrain_sizeState(i) = 0_pInt
do j = 1,maxval(homogenization_Noutput)
select case(homogenization_isostrain_output(j,i))
case('ngrains')
homogenization_isostrain_sizePostResults(i) = &
homogenization_isostrain_sizePostResults(i) + 1
end select
enddo
enddo
return
end subroutine
!*********************************************************************
!* initial homogenization state *
!*********************************************************************
function homogenization_isostrain_stateInit(myInstance)
use prec, only: pReal,pInt
implicit none
!* Definition of variables
integer(pInt), intent(in) :: myInstance
real(pReal), dimension(1) :: homogenization_isostrain_stateInit
homogenization_isostrain_stateInit = 0.0_pReal
return
end function
!********************************************************************
! partition material point def grad onto constituents
!********************************************************************
subroutine homogenization_isostrain_partitionDeformation(&
F, & ! partioned def grad per grain
!
F0, & ! initial partioned def grad per grain
avgF, & ! my average def grad
state, & ! my state
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 (3,3,homogenization_maxNgrains), intent(out) :: F
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: F0
real(pReal), dimension (3,3), intent(in) :: avgF
type(p_vec), intent(in) :: state
integer(pInt), intent(in) :: ip,el
integer(pInt) homID, i
! homID = homogenization_typeInstance(mesh_element(3,el))
forall (i = 1:homogenization_Ngrains(mesh_element(3,el))) &
F(:,:,i) = avgF
return
end subroutine
!********************************************************************
! update the internal state of the homogenization scheme
! and tell whether "done" and "happy" with result
!********************************************************************
function homogenization_isostrain_updateState(&
state, & ! my state
!
P, & ! array of current grain stresses
dPdF, & ! array of current grain stiffnesses
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
implicit none
!* Definition of variables
type(p_vec), intent(inout) :: state
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: P
real(pReal), dimension (3,3,3,3,homogenization_maxNgrains), intent(in) :: dPdF
integer(pInt), intent(in) :: ip,el
! integer(pInt) homID
logical, dimension(2) :: homogenization_isostrain_updateState
! homID = homogenization_typeInstance(mesh_element(3,el))
homogenization_isostrain_updateState = .true. ! homogenization at material point converged (done and happy)
return
end function
!********************************************************************
! derive average stress and stiffness from constituent quantities
!********************************************************************
subroutine homogenization_isostrain_averageStressAndItsTangent(&
avgP, & ! average stress at material point
dAvgPdAvgF, & ! average stiffness at material point
!
P, & ! array of current grain stresses
dPdF, & ! array of current grain stiffnesses
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 (3,3), intent(out) :: avgP
real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: P
real(pReal), dimension (3,3,3,3,homogenization_maxNgrains), intent(in) :: dPdF
integer(pInt), intent(in) :: ip,el
logical homogenization_isostrain_stateUpdate
integer(pInt) homID, i, Ngrains
! homID = homogenization_typeInstance(mesh_element(3,el))
Ngrains = homogenization_Ngrains(mesh_element(3,el))
avgP = sum(P,3)/Ngrains
dAvgPdAvgF = sum(dPdF,5)/Ngrains
return
end subroutine
!********************************************************************
! return array of homogenization results for post file inclusion
!********************************************************************
pure function homogenization_isostrain_postResults(&
state, & ! my state
ip, & ! my integration point
el & ! my element
)
use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_element
use material, only: homogenization_typeInstance,homogenization_Noutput
implicit none
!* Definition of variables
type(p_vec), intent(in) :: state
integer(pInt), intent(in) :: ip,el
integer(pInt) homID,o,c
real(pReal), dimension(homogenization_isostrain_sizePostResults(homogenization_typeInstance(mesh_element(3,el)))) :: &
homogenization_isostrain_postResults
homID = homogenization_typeInstance(mesh_element(3,el))
c = 0_pInt
homogenization_isostrain_postResults = 0.0_pReal
do o = 1,homogenization_Noutput(mesh_element(3,el))
select case(homogenization_isostrain_output(o,homID))
case ('ngrains')
homogenization_isostrain_postResults(c+1) = homogenization_isostrain_Ngrains(homID)
c = c + 1
end select
enddo
return
end function
END MODULE

3
trunk/lattice.f90
View File

@ -566,7 +566,6 @@ subroutine lattice_init()
allocate(lattice_interactionSlipTwin(lattice_maxNslip,lattice_maxNtwin,lattice_Nstructure)); lattice_interactionSlipTwin = 0_pInt allocate(lattice_interactionSlipTwin(lattice_maxNslip,lattice_maxNtwin,lattice_Nstructure)); lattice_interactionSlipTwin = 0_pInt
allocate(lattice_interactionTwinTwin(lattice_maxNtwin,lattice_maxNtwin,lattice_Nstructure)); lattice_interactionTwinTwin = 0_pInt allocate(lattice_interactionTwinTwin(lattice_maxNtwin,lattice_maxNtwin,lattice_Nstructure)); lattice_interactionTwinTwin = 0_pInt
write(6,*) 'lattice Nstructure',lattice_Nstructure
end subroutine end subroutine
@ -595,7 +594,6 @@ function lattice_initializeStructure(struct,CoverA)
integer(pInt) lattice_initializeStructure integer(pInt) lattice_initializeStructure
write(6,*) 'initialize structure', struct
select case(struct(1:3)) ! check first three chars of structure name select case(struct(1:3)) ! check first three chars of structure name
case ('fcc') case ('fcc')
myStructure = 1_pInt myStructure = 1_pInt
@ -705,7 +703,6 @@ function lattice_initializeStructure(struct,CoverA)
endif endif
lattice_initializeStructure = myStructure lattice_initializeStructure = myStructure
write(6,*) 'lattice_initializeStructure', myStructure
end function end function

48
trunk/material.f90
View File

@ -51,8 +51,8 @@ integer(pInt), dimension(:), allocatable :: homogenization_Ngrains, &
integer(pInt), dimension(:,:), allocatable :: microstructure_phase, & ! phase IDs of each microstructure integer(pInt), dimension(:,:), allocatable :: microstructure_phase, & ! phase IDs of each microstructure
microstructure_texture ! texture IDs of each microstructure microstructure_texture ! texture IDs of each microstructure
real(pReal), dimension(:,:), allocatable :: microstructure_fraction ! vol fraction of each constituent in microstructure real(pReal), dimension(:,:), allocatable :: microstructure_fraction ! vol fraction of each constituent in microstructure
real(pReal), dimension(:,:,:), allocatable :: material_volFrac ! vol fraction of grain within phase (?) integer(pInt), dimension(:,:,:), allocatable :: material_volFrac, & ! vol fraction of grain within phase (?)
integer(pInt), dimension(:,:,:), allocatable :: material_phase ! phase of each grain,IP,element material_phase ! phase of each grain,IP,element
real(pReal), dimension(:,:,:,:), allocatable :: material_EulerAngles ! initial orientation of each grain,IP,element real(pReal), dimension(:,:,:,:), allocatable :: material_EulerAngles ! initial orientation of each grain,IP,element
real(pReal), dimension(:,:,:), allocatable :: texture_Gauss, & ! data of each Gauss component real(pReal), dimension(:,:,:), allocatable :: texture_Gauss, & ! data of each Gauss component
texture_Fiber ! data of each Fiber component texture_Fiber ! data of each Fiber component
@ -74,13 +74,9 @@ subroutine material_init()
integer(pInt) i integer(pInt) i
if(.not. IO_open_file(fileunit,material_configFile)) call IO_error (100) ! corrupt config file if(.not. IO_open_file(fileunit,material_configFile)) call IO_error (100) ! corrupt config file
write(6,*) 'parsing homogenization...'
call material_parseHomogenization(fileunit,material_partHomogenization) call material_parseHomogenization(fileunit,material_partHomogenization)
write(6,*) 'parsing microstrcuture...'
call material_parseMicrostructure(fileunit,material_partMicrostructure) call material_parseMicrostructure(fileunit,material_partMicrostructure)
write(6,*) 'parsing texture...'
call material_parseTexture(fileunit,material_partTexture) call material_parseTexture(fileunit,material_partTexture)
write(6,*) 'parsing phase...'
call material_parsePhase(fileunit,material_partPhase) call material_parsePhase(fileunit,material_partPhase)
close(fileunit) close(fileunit)
@ -104,9 +100,7 @@ subroutine material_init()
write (6,'(a32,x,i4)') microstructure_name(i),microstructure_Nconstituents(i) write (6,'(a32,x,i4)') microstructure_name(i),microstructure_Nconstituents(i)
enddo enddo
write(6,*) 'populating grains...'
call material_populateGrains() call material_populateGrains()
write(6,*) 'populating grains finished...'
end subroutine end subroutine
@ -129,12 +123,17 @@ subroutine material_parseHomogenization(file,myPart)
Nsections = IO_countSections(file,myPart) Nsections = IO_countSections(file,myPart)
material_Nhomogenization = Nsections material_Nhomogenization = Nsections
write (6,*) 'homogenization sections found',material_Nhomogenization
allocate(homogenization_name(Nsections)); homogenization_name = '' allocate(homogenization_name(Nsections)); homogenization_name = ''
allocate(homogenization_type(Nsections)); homogenization_type = '' allocate(homogenization_type(Nsections)); homogenization_type = ''
allocate(homogenization_typeInstance(Nsections)); homogenization_typeInstance = 0_pInt allocate(homogenization_typeInstance(Nsections)); homogenization_typeInstance = 0_pInt
allocate(homogenization_Ngrains(Nsections)); homogenization_Ngrains = 0_pInt allocate(homogenization_Ngrains(Nsections)); homogenization_Ngrains = 0_pInt
allocate(homogenization_Noutput(Nsections)); homogenization_Noutput = 0_pInt
write(6,*) 'scanning for (output)',homogenization_Noutput
homogenization_Noutput = IO_countTagInPart(file,myPart,'(output)',Nsections) homogenization_Noutput = IO_countTagInPart(file,myPart,'(output)',Nsections)
write(6,*) 'count of (output)',homogenization_Noutput
rewind(file) rewind(file)
line = '' line = ''
@ -434,6 +433,7 @@ subroutine material_populateGrains()
integer(pInt), dimension (:,:), allocatable :: Ngrains integer(pInt), dimension (:,:), allocatable :: Ngrains
integer(pInt), dimension (microstructure_maxNconstituents) :: NgrainsOfConstituent integer(pInt), dimension (microstructure_maxNconstituents) :: NgrainsOfConstituent
real(pReal), dimension (:,:), allocatable :: volume
real(pReal), dimension (:), allocatable :: volFracOfGrain, phaseOfGrain real(pReal), dimension (:), allocatable :: volFracOfGrain, phaseOfGrain
real(pReal), dimension (:,:), allocatable :: orientationOfGrain real(pReal), dimension (:,:), allocatable :: orientationOfGrain
real(pReal), dimension (3) :: orientation real(pReal), dimension (3) :: orientation
@ -448,8 +448,9 @@ subroutine material_populateGrains()
allocate(material_EulerAngles(3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; material_EulerAngles = 0.0_pReal allocate(material_EulerAngles(3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; material_EulerAngles = 0.0_pReal
allocate(Ngrains(material_Nhomogenization,material_Nmicrostructure)); Ngrains = 0_pInt allocate(Ngrains(material_Nhomogenization,material_Nmicrostructure)); Ngrains = 0_pInt
allocate(volume(material_Nhomogenization,material_Nmicrostructure)); volume = 0.0_pReal
! count grains per homog/micro pair ! count grains and total volume per homog/micro pair
do e = 1,mesh_NcpElems do e = 1,mesh_NcpElems
homog = mesh_element(3,e) homog = mesh_element(3,e)
micro = mesh_element(4,e) micro = mesh_element(4,e)
@ -458,6 +459,7 @@ subroutine material_populateGrains()
if (micro < 1 .or. micro > material_Nmicrostructure) & ! out of bounds if (micro < 1 .or. micro > material_Nmicrostructure) & ! out of bounds
call IO_error(140,e,0,0) call IO_error(140,e,0,0)
Ngrains(homog,micro) = Ngrains(homog,micro) + homogenization_Ngrains(homog) * FE_Nips(mesh_element(2,e)) Ngrains(homog,micro) = Ngrains(homog,micro) + homogenization_Ngrains(homog) * FE_Nips(mesh_element(2,e))
volume(homog,micro) = volume(homog,micro) + sum(mesh_ipVolume(:,e))
enddo enddo
allocate(volFracOfGrain(maxval(Ngrains))) ! reserve memory for maximum case allocate(volFracOfGrain(maxval(Ngrains))) ! reserve memory for maximum case
@ -480,16 +482,15 @@ subroutine material_populateGrains()
do e = 1,mesh_NcpElems ! check each element do e = 1,mesh_NcpElems ! check each element
if (mesh_element(3,e) == homog .and. mesh_element(4,e) == micro) then ! my combination of homog and micro if (mesh_element(3,e) == homog .and. mesh_element(4,e) == micro) then ! my combination of homog and micro
forall (i = 1:FE_Nips(mesh_element(2,e))) & ! loop over IPs forall (i = 1:FE_Nips(mesh_element(2,e))) & ! loop over IPs
volFracOfGrain(grain+(i-1)*dGrains+1:grain+i*dGrains) = mesh_ipVolume(i,e)/dGrains ! assign IPvolfrac/Ngrains to grains volFracOfGrain(grain+(i-1)*dGrains+1:grain+i*dGrains) = &
mesh_ipVolume(i,e)/volume(homog,micro)/dGrains ! assign IPvolfrac/Ngrains to grains
grain = grain + FE_Nips(mesh_element(2,e)) * dGrains ! wind forward by Nips*NgrainsPerIP grain = grain + FE_Nips(mesh_element(2,e)) * dGrains ! wind forward by Nips*NgrainsPerIP
endif endif
enddo enddo
write (6,*) 'now at grain count',grain
! ---------------------------------------------------------------------------- ! ----------------------------------------------------------------------------
NgrainsOfConstituent = 0_pInt NgrainsOfConstituent = 0_pInt
forall (i = 1:microstructure_Nconstituents(micro)) & forall (i = 1:microstructure_Nconstituents(micro)) &
NgrainsOfConstituent(i) = nint(microstructure_fraction(i,micro) * myNgrains, pInt) NgrainsOfConstituent(i) = nint(microstructure_fraction(i,micro) * myNgrains, pInt)
write (6,*) 'NgrainsOfConstituent',NgrainsOfConstituent
do while (sum(NgrainsOfConstituent) /= myNgrains) ! total grain count over constituents wrong? do while (sum(NgrainsOfConstituent) /= myNgrains) ! total grain count over constituents wrong?
sgn = sign(1_pInt, myNgrains - sum(NgrainsOfConstituent)) ! direction of required change sgn = sign(1_pInt, myNgrains - sum(NgrainsOfConstituent)) ! direction of required change
extreme = 0.0_pReal extreme = 0.0_pReal
@ -502,8 +503,6 @@ subroutine material_populateGrains()
enddo enddo
NgrainsOfConstituent(t) = NgrainsOfConstituent(t) + sgn ! change that by one NgrainsOfConstituent(t) = NgrainsOfConstituent(t) + sgn ! change that by one
end do end do
write (6,*) 'fixed NgrainsOfConstituent',NgrainsOfConstituent
! ---------------------------------------------------------------------------- ! ----------------------------------------------------------------------------
phaseOfGrain = 0_pInt phaseOfGrain = 0_pInt
orientationOfGrain = 0.0_pReal orientationOfGrain = 0.0_pReal
@ -515,26 +514,21 @@ subroutine material_populateGrains()
phaseOfGrain(grain+1:grain+NgrainsOfConstituent(i)) = phaseID ! assign resp. phase phaseOfGrain(grain+1:grain+NgrainsOfConstituent(i)) = phaseID ! assign resp. phase
myNorientations = ceiling(float(NgrainsOfConstituent(i))/texture_symmetry(textureID)) ! max number of unique orientations (excl. symmetry) myNorientations = ceiling(float(NgrainsOfConstituent(i))/texture_symmetry(textureID)) ! max number of unique orientations (excl. symmetry)
write (6,'(a32,x,i6,x,f5.3,x,i6)') &
phase_name(phaseID),NgrainsOfConstituent(i),real(NgrainsOfConstituent(i)/myNgrains),myNorientations
constituentGrain = 0_pInt ! constituent grain index constituentGrain = 0_pInt ! constituent grain index
! --------- ! ---------
if (texture_ODFfile(textureID) == '') then ! dealing with texture components if (texture_ODFfile(textureID) == '') then ! dealing with texture components
! --------- ! ---------
do t = 1,texture_Ngauss(textureID) ! loop over Gauss components do t = 1,texture_Ngauss(textureID) ! loop over Gauss components
write (6,*) 'gauss',t,int(myNorientations*texture_Gauss(5,t,textureID))
do g = 1,int(myNorientations*texture_Gauss(5,t,textureID)) ! loop over required grain count do g = 1,int(myNorientations*texture_Gauss(5,t,textureID)) ! loop over required grain count
orientationOfGrain(:,grain+constituentGrain+g) = & orientationOfGrain(:,grain+constituentGrain+g) = &
math_sampleGaussOri(texture_Gauss(1:3,t,textureID),& math_sampleGaussOri(texture_Gauss(1:3,t,textureID),&
texture_Gauss( 4,t,textureID)) texture_Gauss( 4,t,textureID))
enddo enddo
constituentGrain = constituentGrain + int(myNorientations*texture_Gauss(5,t,textureID)) constituentGrain = constituentGrain + int(myNorientations*texture_Gauss(5,t,textureID))
write (6,*) 'now at constituent grain',constituentGrain
enddo enddo
do t = 1,texture_Nfiber(textureID) ! loop over fiber components do t = 1,texture_Nfiber(textureID) ! loop over fiber components
write (6,*) 'fiber',t,int(myNorientations*texture_Fiber(6,t,textureID))
do g = 1,int(myNorientations*texture_Fiber(6,t,textureID)) ! loop over required grain count do g = 1,int(myNorientations*texture_Fiber(6,t,textureID)) ! loop over required grain count
orientationOfGrain(:,grain+constituentGrain+g) = & orientationOfGrain(:,grain+constituentGrain+g) = &
math_sampleFiberOri(texture_Fiber(1:2,t,textureID),& math_sampleFiberOri(texture_Fiber(1:2,t,textureID),&
@ -542,14 +536,11 @@ subroutine material_populateGrains()
texture_Fiber( 5,t,textureID)) texture_Fiber( 5,t,textureID))
enddo enddo
constituentGrain = constituentGrain + int(myNorientations*texture_fiber(6,t,textureID)) constituentGrain = constituentGrain + int(myNorientations*texture_fiber(6,t,textureID))
write (6,*) 'now at constituent grain',constituentGrain
enddo enddo
write (6,*) 'looping',constituentGrain+1,myNorientations
do j = constituentGrain+1,myNorientations ! fill remainder with random do j = constituentGrain+1,myNorientations ! fill remainder with random
orientationOfGrain(:,grain+j) = math_sampleRandomOri() orientationOfGrain(:,grain+j) = math_sampleRandomOri()
enddo enddo
write (6,*) 'done...'
! --------- ! ---------
else ! hybrid IA else ! hybrid IA
! --------- ! ---------
@ -592,9 +583,6 @@ subroutine material_populateGrains()
!exchange in MC steps to improve result... !exchange in MC steps to improve result...
! ---------------------------------------------------------------------------- ! ----------------------------------------------------------------------------
!write(6,*) ''
!write(6,*) 'USER DEFINED OUTPUT'
!write(6,'(7(a10,x),a10)') 'element','ip','Ngrains','volFrac','phase','phi1','Phi','phi2'
grain = 0_pInt ! microstructure grain index grain = 0_pInt ! microstructure grain index
do e = 1,mesh_NcpElems ! check each element do e = 1,mesh_NcpElems ! check each element
if (mesh_element(3,e) == homog .and. mesh_element(4,e) == micro) then ! my combination of homog and micro if (mesh_element(3,e) == homog .and. mesh_element(4,e) == micro) then ! my combination of homog and micro
@ -603,13 +591,9 @@ subroutine material_populateGrains()
material_phase(g,i,e) = phaseOfGrain(grain+(i-1)*dGrains+g) material_phase(g,i,e) = phaseOfGrain(grain+(i-1)*dGrains+g)
material_EulerAngles(:,g,i,e) = orientationOfGrain(:,grain+(i-1)*dGrains+g) material_EulerAngles(:,g,i,e) = orientationOfGrain(:,grain+(i-1)*dGrains+g)
end forall end forall
!do i = 1,FE_Nips(mesh_element(2,e)) write (6,*) e
! write(6,'(3(i10,x),e10.3,x,i10,x,3(f10.1,x))') e, i, dGrains, sum(material_volFrac(:,i,e)), & write (6,*) material_phase(:,:,e)
! sum(material_phase(:,i,e)), & write (6,*) material_EulerAngles(:,:,:,e)
! sum(material_EulerAngles(1,:,i,e)), &
! sum(material_EulerAngles(2,:,i,e)), &
! sum(material_EulerAngles(3,:,i,e))
!end do
grain = grain + FE_Nips(mesh_element(2,e)) * dGrains ! wind forward by Nips*NgrainsPerIP grain = grain + FE_Nips(mesh_element(2,e)) * dGrains ! wind forward by Nips*NgrainsPerIP
endif endif
enddo enddo

7
trunk/math.f90
View File

@ -1522,15 +1522,12 @@ math_sampleFiberOri = math_RtoEuler(math_mul33x33(pRot,math_mul33x33(fRot,oRot))
real(pReal) FE(3,3),R(3,3),U(3,3),CE(3,3),EW1,EW2,EW3,EB1(3,3),EB2(3,3),EB3(3,3),UI(3,3),det real(pReal) FE(3,3),R(3,3),U(3,3),CE(3,3),EW1,EW2,EW3,EB1(3,3),EB2(3,3),EB3(3,3),UI(3,3),det
error = .false. error = .false.
!!$OMP CRITICAL (evilmatmul) ce = math_mul33x33(transpose(FE),FE)
ce=math_mul33x33(transpose(fe),fe)
!!$OMP END CRITICAL (evilmatmul)
CALL math_spectral1(CE,EW1,EW2,EW3,EB1,EB2,EB3) CALL math_spectral1(CE,EW1,EW2,EW3,EB1,EB2,EB3)
U=DSQRT(EW1)*EB1+DSQRT(EW2)*EB2+DSQRT(EW3)*EB3 U=DSQRT(EW1)*EB1+DSQRT(EW2)*EB2+DSQRT(EW3)*EB3
call math_invert3x3(U,UI,det,error) call math_invert3x3(U,UI,det,error)
if (.not. error) R = math_mul33x33(fe,ui) if (.not. error) R = math_mul33x33(FE,UI)
return return

37
trunk/mesh.f90
View File

@ -1009,7 +1009,6 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
END SUBROUTINE END SUBROUTINE
!******************************************************************** !********************************************************************
! get count of elements, nodes, and cp elements in mesh ! get count of elements, nodes, and cp elements in mesh
! for subsequent array allocations ! for subsequent array allocations
@ -1685,18 +1684,18 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit)
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
!write (6,*) write (6,*)
!write (6,*) "Input Parser: IP NEIGHBORHOOD" write (6,*) "Input Parser: IP NEIGHBORHOOD"
!write (6,*) write (6,*)
!write (6,"(a10,x,a10,x,a10,x,a3,x,a13,x,a13)") "elem","IP","neighbor","","elemNeighbor","ipNeighbor" write (6,"(a10,x,a10,x,a10,x,a3,x,a13,x,a13)") "elem","IP","neighbor","","elemNeighbor","ipNeighbor"
!do e = 1,mesh_NcpElems ! loop over cpElems do e = 1,mesh_NcpElems ! loop over cpElems
! t = mesh_element(2,e) ! get elemType t = mesh_element(2,e) ! get elemType
! do i = 1,FE_Nips(t) ! loop over IPs of elem do i = 1,FE_Nips(t) ! loop over IPs of elem
! do n = 1,FE_NipNeighbors(t) ! loop over neighbors of IP do n = 1,FE_NipNeighbors(t) ! loop over neighbors of IP
! write (6,"(i10,x,i10,x,i10,x,a3,x,i13,x,i13)") e,i,n,'-->',mesh_ipNeighborhood(1,n,i,e),mesh_ipNeighborhood(2,n,i,e) write (6,"(i10,x,i10,x,i10,x,a3,x,i13,x,i13)") e,i,n,'-->',mesh_ipNeighborhood(1,n,i,e),mesh_ipNeighborhood(2,n,i,e)
! enddo enddo
! enddo enddo
!enddo enddo
write (6,*) write (6,*)
write (6,*) "Input Parser: ELEMENT VOLUME" write (6,*) "Input Parser: ELEMENT VOLUME"
write (6,*) write (6,*)
@ -1706,9 +1705,9 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit)
do e = 1,mesh_NcpElems do e = 1,mesh_NcpElems
do i = 1,FE_Nips(mesh_element(2,e)) do i = 1,FE_Nips(mesh_element(2,e))
write (6,"(i5,x,i5,x,e15.8)") e,i,mesh_IPvolume(i,e) write (6,"(i5,x,i5,x,e15.8)") e,i,mesh_IPvolume(i,e)
! do f = 1,FE_NipNeighbors(mesh_element(2,e)) do f = 1,FE_NipNeighbors(mesh_element(2,e))
! write (6,"(i33,x,e15.8,x,3(f6.3,x))") f,mesh_ipArea(f,i,e),mesh_ipAreaNormal(:,f,i,e) write (6,"(i33,x,e15.8,x,3(f6.3,x))") f,mesh_ipArea(f,i,e),mesh_ipAreaNormal(:,f,i,e)
! end do end do
end do end do
end do end do
!write (6,*) !write (6,*)
@ -1746,9 +1745,9 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit)
write (6,*) mesh_maxValStateVar(1), " : maximum homogenization index" write (6,*) mesh_maxValStateVar(1), " : maximum homogenization index"
write (6,*) mesh_maxValStateVar(2), " : maximum microstructure index" write (6,*) mesh_maxValStateVar(2), " : maximum microstructure index"
write (6,*) write (6,*)
write (fmt,"(a,i5,a)") "(9(x),a1,x,",mesh_maxValStateVar(2),"(i8))" write (fmt,"(a,i5,a)") "(9(x),a2,x,",mesh_maxValStateVar(2),"(i8))"
write (6,fmt) "+",math_range(mesh_maxValStateVar(2)) write (6,fmt) "+-",math_range(mesh_maxValStateVar(2))
write (fmt,"(a,i5,a)") "(i8,x,a1,x,",mesh_maxValStateVar(2),"(i8))" write (fmt,"(a,i5,a)") "(i8,x,a2,x,",mesh_maxValStateVar(2),"(i8))"
do i=1,mesh_maxValStateVar(1) ! loop over all (possibly assigned) homogenizations do i=1,mesh_maxValStateVar(1) ! loop over all (possibly assigned) homogenizations
write (6,fmt) i,"| ",mesh_HomogMicro(i,:) ! loop over all (possibly assigned) microstrcutures write (6,fmt) i,"| ",mesh_HomogMicro(i,:) ! loop over all (possibly assigned) microstrcutures
enddo enddo

160
trunk/mpie_cpfem_marc.f90
View File

@ -37,19 +37,18 @@
include "math.f90" ! uses prec include "math.f90" ! uses prec
include "IO.f90" ! uses prec, debug, math include "IO.f90" ! uses prec, debug, math
include "FEsolving.f90" ! uses prec, IO include "FEsolving.f90" ! uses prec, IO
include "mesh.f90" ! uses prec, IO, math, FEsolving include "mesh.f90" ! uses prec, math, IO, FEsolving
include "material.f90" ! uses prec, math, IO, mesh include "material.f90" ! uses prec, math, IO, mesh
include "lattice.f90" ! uses prec, math include "lattice.f90" ! uses prec, math, IO, material
include "constitutive_phenomenological.f90" ! uses prec, math, IO, lattice, material, debug include "constitutive_phenomenological.f90" ! uses prec, math, IO, lattice, material, debug
include "constitutive_j2.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_dislobased.f90" ! uses prec, math, IO, lattice, material, debug
include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug
include "crystallite.f90" ! uses
include "homogenization_isostrain.f90" ! uses
include "homogenization.f90" ! uses
include "CPFEM.f90" ! uses prec, math, mesh, constitutive, FEsolving, debug, lattice, IO, crystallite include "CPFEM.f90" ! uses prec, math, mesh, constitutive, FEsolving, debug, lattice, IO, crystallite
SUBROUTINE hypela2(d,g,e,de,s,t,dt,ngens,n,nn,kcus,matus,ndi,&
nshear,disp,dispt,coord,ffn,frotn,strechn,eigvn,ffn1,&
frotn1,strechn1,eigvn1,ncrd,itel,ndeg,ndm,&
nnode,jtype,lclass,ifr,ifu)
!******************************************************************** !********************************************************************
! This is the Marc material routine ! This is the Marc material routine
@ -75,51 +74,6 @@
! update+finite+large disp+constant d) in the parameter section of ! update+finite+large disp+constant d) in the parameter section of
! input deck. ! input deck.
! !
!
! d stress strain law to be formed
! g change in stress due to temperature effects
! e total elastic strain
! de increment of strain
! s stress - should be updated by user
! t state variables (comes in at t=n, must be updated
! to have state variables at t=n+1)
! dt increment of state variables
! ngens size of stress - strain law
! n element number
! nn integration point number
! kcus(1) layer number
! kcus(2) internal layer number
! matus(1) user material identification number
! matus(2) internal material identification number
! ndi number of direct components
! nshear number of shear components
! disp incremental displacements
! dispt displacements at t=n (at assembly, lovl=4) and
! displacements at t=n+1 (at stress recovery, lovl=6)
! coord coordinates
! ncrd number of coordinates
! ndeg number of degrees of freedom
! itel dimension of F and R, either 2 or 3
! nnode number of nodes per element
! jtype element type
! lclass element class
! ifr set to 1 if R has been calculated
! ifu set to 1 if strech has been calculated
!
! at t=n :
!
! ffn deformation gradient
! frotn rotation tensor
! strechn square of principal stretch ratios, lambda(i)
! eigvn(i,j) i principal direction components for j eigenvalues
!
! at t=n+1 :
!
! ffn1 deformation gradient
! frotn1 rotation tensor
! strechn1 square of principal stretch ratios, lambda(i)
! eigvn1(i,j) i principal direction components for j eigenvalues
!
! The following operation obtains U (stretch tensor) at t=n+1 : ! The following operation obtains U (stretch tensor) at t=n+1 :
! !
! call scla(un1,0.d0,itel,itel,1) ! call scla(un1,0.d0,itel,itel,1)
@ -131,12 +85,49 @@
!2 continue !2 continue
!3 continue !3 continue
! !
!********************************************************************
subroutine hypela2(&
d,& ! stress strain law to be formed
g,& ! change in stress due to temperature effects
e,& ! total elastic strain
de,& ! increment of strain
s,& ! stress - should be updated by user
t,& ! state variables (comes in at t=n, must be updated to have state variables at t=n+1)
dt,& ! increment of state variables
ngens,& ! size of stress - strain law
n,& ! element number
nn,& ! integration point number
kcus,& ! (1) layer number, (2) internal layer number
matus,& ! (1) user material identification number, (2) internal material identification number
ndi,& ! number of direct components
nshear,& ! number of shear components
disp,& ! incremental displacements
dispt,& ! displacements at t=n (at assembly, lovl=4) and displacements at t=n+1 (at stress recovery, lovl=6)
coord,& ! coordinates
ffn,& ! deformation gradient
frotn,& ! rotation tensor
strechn,& ! square of principal stretch ratios, lambda(i)
eigvn,& ! i principal direction components for j eigenvalues
ffn1,& ! deformation gradient
frotn1,& ! rotation tensor
strechn1,& ! square of principal stretch ratios, lambda(i)
eigvn1,& ! i principal direction components for j eigenvalues
ncrd,& ! number of coordinates
itel,& ! dimension of F and R, either 2 or 3
ndeg,& ! number of degrees of freedom ==> is this at correct list position ?!?
ndm,& !
nnode,& ! number of nodes per element
jtype,& ! element type
lclass,& ! element class
ifr,& ! set to 1 if R has been calculated
ifu & ! set to 1 if stretch has been calculated
)
use prec, only: pReal,pInt, ijaco use prec, only: pReal,pInt, ijaco
use FEsolving use FEsolving
use CPFEM, only: CPFEM_general use CPFEM, only: CPFEM_general
use math, only: invnrmMandel use math, only: invnrmMandel
use debug use debug
!
implicit none implicit none
! ** Start of generated type statements ** ! ** Start of generated type statements **
@ -146,7 +137,7 @@
integer(pInt) ndi, ndm, ngens, nn, nnode, nshear integer(pInt) ndi, ndm, ngens, nn, nnode, nshear
real(pReal) s, strechn, strechn1, t real(pReal) s, strechn, strechn1, t
! ** End of generated type statements ** ! ** End of generated type statements **
!
dimension e(*),de(*),t(*),dt(*),g(*),d(ngens,*),s(*), n(2),coord(ncrd,*),disp(ndeg,*),matus(2),dispt(ndeg,*),ffn(itel,*),& dimension e(*),de(*),t(*),dt(*),g(*),d(ngens,*),s(*), n(2),coord(ncrd,*),disp(ndeg,*),matus(2),dispt(ndeg,*),ffn(itel,*),&
frotn(itel,*),strechn(itel),eigvn(itel,*),ffn1(itel,*),frotn1(itel,*),strechn1(itel),eigvn1(itel,*),kcus(2), lclass(2) frotn(itel,*),strechn(itel),eigvn(itel,*),ffn1(itel,*),frotn1(itel,*),strechn1(itel),eigvn1(itel,*),kcus(2), lclass(2)
@ -158,8 +149,6 @@
integer(pInt) computationMode,i integer(pInt) computationMode,i
! write(6,'(3(3(f10.3,x),/))') ffn1(:,1),ffn1(:,2),ffn1(:,3)
if (inc == 0) then if (inc == 0) then
cycleCounter = 4 cycleCounter = 4
else else
@ -169,13 +158,7 @@
outdatedFFN1 = .false. outdatedFFN1 = .false.
write (6,*) n(1),nn,'cycleCounter',cycleCounter write (6,*) n(1),nn,'cycleCounter',cycleCounter
call debug_info() ! output of debugging/performance statistics of former call debug_info() ! output of debugging/performance statistics of former
debug_cutbackDistribution = 0_pInt ! initialize debugging data call debug_reset()
debug_InnerLoopDistribution = 0_pInt
debug_OuterLoopDistribution = 0_pInt
debug_cumLpTicks = 0
debug_cumDotStateTicks = 0
debug_cumLpCalls = 0_pInt
debug_cumDotStateCalls = 0_pInt
endif endif
endif endif
if (cptim > theTime .or. theInc /= inc) then ! reached convergence if (cptim > theTime .or. theInc /= inc) then ! reached convergence
@ -213,48 +196,47 @@
return return
END SUBROUTINE end subroutine
!
SUBROUTINE plotv(v,s,sp,etot,eplas,ecreep,t,m,nn,layer,ndi,nshear,jpltcd)
!******************************************************************** !********************************************************************
! This routine sets user defined output variables for Marc ! This routine sets user defined output variables for Marc
!******************************************************************** !********************************************************************
! !
! select a variable contour plotting (user subroutine). ! select a variable contour plotting (user subroutine).
! !
! v variable
! s (idss) stress array
! sp stresses in preferred direction
! etot total strain (generalized)
! eplas total plastic strain
! ecreep total creep strain
! t current temperature
! m element number
! nn integration point number
! layer layer number
! ndi (3) number of direct stress components
! nshear (3) number of shear stress components
!
!******************************************************************** !********************************************************************
subroutine plotv(&
v,& ! variable
s,& ! stress array
sp,& ! stresses in preferred direction
etot,& ! total strain (generalized)
eplas,& ! total plastic strain
ecreep,& ! total creep strain
t,& ! current temperature
m,& ! element number
nn,& ! integration point number
layer,& ! layer number
ndi,& ! number of direct stress components
nshear,& ! number of shear stress components
jpltcd & ! user variable index
)
use prec, only: pReal,pInt use prec, only: pReal,pInt
use CPFEM, only: CPFEM_results, CPFEM_Nresults
use constitutive, only: constitutive_maxSizePostResults
use mesh, only: mesh_FEasCP use mesh, only: mesh_FEasCP
use homogenization, only: materialpoint_results
implicit none implicit none
!
real(pReal) s(*),etot(*),eplas(*),ecreep(*),sp(*) real(pReal) s(*),etot(*),eplas(*),ecreep(*),sp(*)
real(pReal) v, t(*) real(pReal) v, t(*)
integer(pInt) m, nn, layer, ndi, nshear, jpltcd integer(pInt) m, nn, layer, ndi, nshear, jpltcd
!
! assign result variable v = materialpoint_results(jpltcd,nn,mesh_FEasCP('elem', m))
v = CPFEM_results(mod(jpltcd-1_pInt, CPFEM_Nresults+constitutive_maxSizePostResults)+1_pInt,&
(jpltcd-1_pInt)/(CPFEM_Nresults+constitutive_maxSizePostResults)+1_pInt,&
nn, mesh_FEasCP('elem', m))
return return
END SUBROUTINE
! end subroutine
!
! subroutine utimestep(timestep,timestepold,icall,time,timeloadcase) ! subroutine utimestep(timestep,timestepold,icall,time,timeloadcase)
!******************************************************************** !********************************************************************
! This routine modifies the addaptive time step of Marc ! This routine modifies the addaptive time step of Marc

15
trunk/prec.f90
View File

@ -20,14 +20,15 @@
! *** Numerical parameters *** ! *** Numerical parameters ***
integer(pInt), parameter :: ijaco = 1_pInt ! frequency of FEM Jacobi update integer(pInt), parameter :: ijaco = 1_pInt ! frequency of FEM Jacobi update
integer(pInt), parameter :: nCutback = 20_pInt ! max cutbacks accounted for in debug distribution
integer(pInt), parameter :: nReg = 1_pInt ! regularization attempts for Jacobi inversion
real(pReal), parameter :: pert_Fg = 1.0e-6_pReal ! strain perturbation for FEM Jacobi real(pReal), parameter :: pert_Fg = 1.0e-6_pReal ! strain perturbation for FEM Jacobi
integer(pInt), parameter :: nOuter = 20_pInt ! outer loop limit 20 integer(pInt), parameter :: nReg = 1_pInt ! regularization attempts for Jacobi inversion
integer(pInt), parameter :: nInner = 200_pInt ! inner loop limit 200 integer(pInt), parameter :: nHomog = 10_pInt ! homogenization loop limit
real(pReal), parameter :: reltol_Outer = 1.0e-5_pReal ! relative tolerance in outer loop (state) integer(pInt), parameter :: nCryst = 10_pInt ! crystallite loop limit (state update)
real(pReal), parameter :: reltol_Inner = 1.0e-6_pReal ! relative tolerance in inner loop (Lp) integer(pInt), parameter :: nStress = 20_pInt ! stress loop limit
real(pReal), parameter :: abstol_Inner = 1.0e-8_pReal ! absolute tolerance in inner loop (Lp) integer(pInt), parameter :: nLp = 40_pInt ! stress loop limit
real(pReal), parameter :: rTol_crystalliteState = 1.0e-5_pReal ! relative tolerance in crystallite state loop
real(pReal), parameter :: rTol_crystalliteStress = 1.0e-6_pReal ! relative tolerance in crystallite stress loop
real(pReal), parameter :: aTol_crystalliteStress = 1.0e-8_pReal ! absolute tolerance in crystallite stress loop
! !
real(pReal), parameter :: resToler = 1.0e-4_pReal ! relative tolerance of residual in GIA iteration real(pReal), parameter :: resToler = 1.0e-4_pReal ! relative tolerance of residual in GIA iteration
real(pReal), parameter :: resAbsol = 1.0e+2_pReal ! absolute tolerance of residual in GIA iteration (corresponds to ~1 Pa) real(pReal), parameter :: resAbsol = 1.0e+2_pReal ! absolute tolerance of residual in GIA iteration (corresponds to ~1 Pa)