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

This is a major update. 

I restructured the subroutine in view of non local constitutive models.
ATTENTION OpenMP parallelization does not work with this version, needs some more work.
ATTENTION CPFEM_GIA8.f90 is not addopted to the new scheme and will NOT work, should actually be renamed CPFEM_RGC.f90 ;-)
I removed crystallite.f90, the routines are now in CPFEM_*.f90.
I removed the marc2005 interface routine since it is outdated.
This commit is contained in:
Franz Roters 2008-11-28 07:39:39 +00:00
parent e8cd7cef94
commit 7b076185b2
9 changed files with 1732 additions and 1176 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1246
trunk/CPFEM_Taylor.f90
View File

File diff suppressed because it is too large Load Diff

902
trunk/CPFEM_Taylor_sequential.f90
View File

@ -9,23 +9,34 @@
! **************************************************************** ! ****************************************************************
! *** 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_Temperature
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_ffn_bar real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_ffn_bar !average FFN per IP
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_ffn1_bar real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_ffn !individual FFN per grain
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_PK1_bar real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_ffn1_bar !average FFN1 per IP
real(pReal), dimension (:,:,:,:,:,:),allocatable :: CPFEM_dPdF_bar real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_ffn1 !individual FFN1 per grain
real(pReal), dimension (:,:,:), allocatable :: CPFEM_stress_bar real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_PK1_bar !average PK1 per IP
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_jaco_bar real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_PK1 !individual PK1 per grain
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_jaco_knownGood real(pReal), dimension (:,:,:,:,:,:), allocatable :: CPFEM_dPdF_bar !average dPdF per IP
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_results real(pReal), dimension (:,:,:,:,:,:), allocatable :: CPFEM_dPdF_bar_old !old average dPdF per IP
real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_Lp real(pReal), dimension (:,:,:,:,:,:,:),allocatable :: CPFEM_dPdF !individual dPdF per grain
real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_Fp_old real(pReal), dimension (:,:,:), allocatable :: CPFEM_stress_bar
real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_Fp_new 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_Fe1
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_Tstar_v
real(pReal), parameter :: CPFEM_odd_stress = 1e15_pReal, CPFEM_odd_jacobian = 1e50_pReal real(pReal), parameter :: CPFEM_odd_stress = 1e15_pReal, CPFEM_odd_jacobian = 1e50_pReal
integer(pInt) :: CPFEM_Nresults = 4_pInt ! three Euler angles plus volume fraction integer(pInt) :: CPFEM_Nresults = 4_pInt ! three Euler angles plus volume fraction
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
logical :: CPFEM_results_aged = .false. ! remember whether results have been aged at inc start logical :: CPFEM_results_aged = .false. ! remember whether results have been aged at inc start
! *** Solution at single crystallite level ***
!
logical, dimension (:,:,:),allocatable :: crystallite_converged !individual covergence flag per grain
! !
CONTAINS CONTAINS
! !
@ -46,12 +57,18 @@
integer(pInt) e,i,g integer(pInt) e,i,g
! !
! *** mpie.marc parameters *** ! *** mpie.marc parameters ***
allocate(CPFEM_Temperature (mesh_maxNips,mesh_NcpElems)) ; CPFEM_Temperature = Temperature allocate(CPFEM_Temperature(mesh_maxNips,mesh_NcpElems)) ; CPFEM_Temperature = Temperature
allocate(CPFEM_ffn_bar (3,3,mesh_maxNips,mesh_NcpElems)) 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 forall(e=1:mesh_NcpElems,i=1:mesh_maxNips) CPFEM_ffn_bar(:,:,i,e) = math_I3
allocate(CPFEM_ffn1_bar (3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_ffn1_bar = CPFEM_ffn_bar allocate(CPFEM_ffn(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems))
allocate(CPFEM_PK1_bar (3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_PK1_bar = 0.0_pReal forall(g=1:constitutive_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,constitutive_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,constitutive_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(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,constitutive_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_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_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 allocate(CPFEM_jaco_knownGood(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_jaco_knownGood = 0.0_pReal
@ -61,31 +78,44 @@
CPFEM_results = 0.0_pReal CPFEM_results = 0.0_pReal
! !
! *** Plastic velocity gradient *** ! *** Plastic velocity gradient ***
allocate(CPFEM_Lp(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Lp = 0.0_pReal allocate(CPFEM_Lp_old(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Lp_old = 0.0_pReal
allocate(CPFEM_Lp_new(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Lp_new = 0.0_pReal
! *** Plastic deformation gradient at (t=t0) and (t=t1) *** ! *** Plastic deformation gradient at (t=t0) and (t=t1) ***
allocate(CPFEM_Fp_new(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Fp_new = 0.0_pReal allocate(CPFEM_Fp_new(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Fp_new = 0.0_pReal
allocate(CPFEM_Fp_old(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) allocate(CPFEM_Fp_old(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems))
forall (e=1:mesh_NcpElems,i=1:mesh_maxNips,g=1:constitutive_maxNgrains) & forall (e=1:mesh_NcpElems,i=1:mesh_maxNips,g=1:constitutive_maxNgrains) &
CPFEM_Fp_old(:,:,g,i,e) = math_EulerToR(constitutive_EulerAngles(:,g,i,e)) ! plastic def gradient reflects init orientation CPFEM_Fp_old(:,:,g,i,e) = math_EulerToR(constitutive_EulerAngles(:,g,i,e)) ! plastic def gradient reflects init orientation
! *** Elastic deformation gradient at (t=t1) ***
allocate(CPFEM_Fe1(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Fe1 = 0.0_pReal
! *** Stress vector at (t=t1) ***
allocate(CPFEM_Tstar_v(6,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Tstar_v = 0.0_pReal
! !
! *** 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 Initialization'
write(6,*) write(6,*)
write(6,*) 'CPFEM_Temperature: ', shape(CPFEM_Temperature) write(6,*) 'CPFEM_Temperature: ', shape(CPFEM_Temperature)
write(6,*) 'CPFEM_ffn_bar: ', shape(CPFEM_ffn_bar) write(6,*) 'CPFEM_ffn_bar: ', shape(CPFEM_ffn_bar)
write(6,*) 'CPFEM_ffn1_bar: ', shape(CPFEM_ffn1_bar) write(6,*) 'CPFEM_ffn: ', shape(CPFEM_ffn)
write(6,*) 'CPFEM_PK1_bar: ', shape(CPFEM_PK1_bar) write(6,*) 'CPFEM_ffn1_bar: ', shape(CPFEM_ffn1_bar)
write(6,*) 'CPFEM_dPdF_bar: ', shape(CPFEM_dPdF_bar) write(6,*) 'CPFEM_ffn1: ', shape(CPFEM_ffn1)
write(6,*) 'CPFEM_stress_bar: ', shape(CPFEM_stress_bar) write(6,*) 'CPFEM_PK1_bar: ', shape(CPFEM_PK1_bar)
write(6,*) 'CPFEM_jaco_bar: ', shape(CPFEM_jaco_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_jaco_knownGood: ', shape(CPFEM_jaco_knownGood)
write(6,*) 'CPFEM_results: ', shape(CPFEM_results) write(6,*) 'CPFEM_results: ', shape(CPFEM_results)
write(6,*) 'CPFEM_Lp: ', shape(CPFEM_Lp) write(6,*) 'CPFEM_Lp_old: ', shape(CPFEM_Lp_old)
write(6,*) 'CPFEM_Fp_old: ', shape(CPFEM_Fp_old) write(6,*) 'CPFEM_Lp_new: ', shape(CPFEM_Lp_new)
write(6,*) 'CPFEM_Fp_new: ', shape(CPFEM_Fp_new) write(6,*) 'CPFEM_Fp_old: ', shape(CPFEM_Fp_old)
write(6,*) 'CPFEM_Fp_new: ', shape(CPFEM_Fp_new)
write(6,*) 'CPFEM_Fe1: ', shape(CPFEM_Fe1)
write(6,*) 'CPFEM_Tstar_v: ', shape(CPFEM_Tstar_v)
write(6,*) write(6,*)
call flush(6) call flush(6)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
@ -142,6 +172,7 @@
call mesh_init() call mesh_init()
call lattice_init() call lattice_init()
call constitutive_init() call constitutive_init()
call crystallite_init()
call CPFEM_init(Temperature) call CPFEM_init(Temperature)
CPFEM_init_done = .true. CPFEM_init_done = .true.
endif endif
@ -186,9 +217,9 @@
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)*CPFEM_ffn1_bar(i,m,CPFEM_in,cp_en)*CPFEM_PK1_bar(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_bar(j,l) + math_I3(j,l)*Kirchhoff_bar(i,k) + &
math_I3(i,l)*Kirchhoff_bar(j,k) + math_I3(j,k)*Kirchhoff_bar(i,l)) math_I3(i,l)*Kirchhoff_bar(j,k) + math_I3(j,k)*Kirchhoff_bar(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_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))
CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel3333to66(J_inverse*H_bar) CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel3333to66(J_inverse*H_bar)
case (3) ! *** collect and return odd result *** case (3) ! *** collect and return odd result ***
CPFEM_Temperature(CPFEM_in,cp_en) = Temperature CPFEM_Temperature(CPFEM_in,cp_en) = Temperature
@ -217,111 +248,704 @@
END SUBROUTINE END SUBROUTINE
! !
! !
!********************************************************** !**********************************************************
!*** calculate the material point behaviour *** !*** calculate the material point behaviour ***
!********************************************************** !**********************************************************
SUBROUTINE CPFEM_MaterialPoint(& SUBROUTINE CPFEM_MaterialPoint(&
updateJaco,& ! flag to initiate Jacobian updating updateJaco,& ! flag to initiate Jacobian updating
CPFEM_dt,& ! Time increment (dt) CPFEM_dt,& ! Time increment (dt)
CPFEM_in,& ! Integration point number CPFEM_in,& ! Integration point number
cp_en) ! Element number cp_en) ! Element number
! !
use prec use prec
use FEsolving, only: theCycle use FEsolving, only: theCycle
use debug use debug
use math, only: math_pDecomposition,math_RtoEuler,inDeg,math_I3,math_invert3x3,math_permut,math_invert,math_delta use math, only: math_pDecomposition,math_RtoEuler,inDeg,math_I3,math_invert3x3,math_permut,math_invert,math_delta
use IO, only: IO_error use IO, only: IO_error
use mesh, only: mesh_element use mesh, only: mesh_element, mesh_NcpElems, FE_Nips
use crystallite ! use crystallite
use constitutive use constitutive
implicit none implicit none
! !
character(len=128) msg integer(pInt) cp_en,CPFEM_in,g,i,e
integer(pInt) cp_en,CPFEM_in,grain integer(pInt) el_start, el_end, ip_start, ip_end
logical updateJaco,error logical updateJaco,error
real(pReal) CPFEM_dt,volfrac real(pReal) CPFEM_dt,volfrac
real(pReal), dimension(3,3) :: U,R,Fe1 real(pReal), dimension(3,3) :: U,R !,Fe1
real(pReal), dimension(3,3) :: PK1 ! real(pReal), dimension(3,3) :: PK1
real(pReal), dimension(3,3,3,3) :: dPdF,dPdF_bar_old ! real(pReal), dimension(3,3,3,3) :: dPdF,dPdF_bar_old
! !
CPFEM_PK1_bar(:,:,CPFEM_in,cp_en) = 0.0_pReal ! zero out average first PK stress CPFEM_PK1_bar = 0.0_pReal ! zero out average first PK stress
if (updateJaco) then !initialize element loop
dPdF_bar_old = CPFEM_dPdF_bar(:,:,:,:,CPFEM_in,cp_en) ! remember former average consistent tangent if (cp_en /= 0_pInt) then
CPFEM_dPdF_bar(:,:,:,:,CPFEM_in,cp_en) = 0.0_pReal ! zero out avg consistent tangent for later assembly el_start = cp_en
endif el_end = cp_en
else
do grain = 1,texture_Ngrains(mesh_element(4,cp_en)) el_start = 1_pInt
dPdF = dPdF_bar_old ! preguess consistent tangent of grain with avg el_end = mesh_NcpElems
endif
if (debugger) then ! prescribe FFN and FFN1 depending on homogenization scheme
!$OMP CRITICAL (write2out) !$OMP PARALLEL DO
write (6,*) 'single crystallite integrating.',cp_en,CPFEM_in,grain do e=el_start,el_end
write (6,'(a,/,3(3(f12.7,x)/))') 'Fg',CPFEM_ffn1_bar(:,:,CPFEM_in,cp_en) if(CPFEM_in /= 0_pInt) then
write (6,'(a,/,3(3(f12.7,x)/))') 'Lp (guess)',CPFEM_Lp(1:3,:,grain,CPFEM_in,cp_en) ip_start = CPFEM_in
write (6,'(a,/,3(3(f12.7,x)/))') 'Fp (old)',CPFEM_Fp_old(1:3,:,grain,CPFEM_in,cp_en) ip_end = CPFEM_in
write (6,'(a,/,3(4(f9.3,x)/))') 'state (old) / MPa',constitutive_state_old(:,grain,CPFEM_in,cp_en)/1e6_pReal else
write (6,'(a,/,3(4(f9.3,x)/))') 'state (new) / MPa',constitutive_state_new(:,grain,CPFEM_in,cp_en)/1e6_pReal ip_start = 1
write (6,*) ip_end = FE_Nips(mesh_element(2,e)) ! iterate over all IPs of this element's type
!$OMP END CRITICAL (write2out) endif
endif do i=ip_start,ip_end
do g=1,texture_Ngrains(mesh_element(4,e))
call SingleCrystallite(msg,PK1,dPdF,& CPFEM_ffn(:,:,g,i,e) = CPFEM_ffn_bar(:,:,i,e) !Taylor homogenization
CPFEM_results(CPFEM_Nresults+1:CPFEM_Nresults+constitutive_Nresults(grain,CPFEM_in,cp_en),& CPFEM_ffn1(:,:,g,i,e) = CPFEM_ffn1_bar(:,:,i,e) !Taylor homogenization
grain,CPFEM_in,cp_en),& end do
CPFEM_Lp(:,:,grain,CPFEM_in,cp_en),& end do
CPFEM_Fp_new(:,:,grain,CPFEM_in,cp_en),Fe1,constitutive_state_new(:,grain,CPFEM_in,cp_en),& ! output up to here end do
CPFEM_dt,cp_en,CPFEM_in,grain,updateJaco,& !$OMP END PARALLEL DO
CPFEM_Temperature(CPFEM_in,cp_en),& ! calculate stress, update state and update jacobian in case needed for all or one ip
CPFEM_ffn1_bar(:,:,CPFEM_in,cp_en),CPFEM_ffn_bar(:,:,CPFEM_in,cp_en),& if (updateJaco) then
CPFEM_Fp_old(:,:,grain,CPFEM_in,cp_en),constitutive_state_old(:,grain,CPFEM_in,cp_en)) CPFEM_dPdF_bar_old = CPFEM_dPdF_bar ! remember former average consistent tangent
CPFEM_dPdF_bar = 0.0_pReal ! zero out avg consistent tangent for later assembly
if (msg /= 'ok') then ! solution not reached --> exit endif
!$OMP CRITICAL (write2out) call SingleCrystallite(updateJaco,CPFEM_dt,el_start,el_end,CPFEM_in)
write(6,*) 'grain loop failed to converge @ EL:',cp_en,' IP:',CPFEM_in !******************************************************************************************************
!$OMP END CRITICAL (write2out) ! check convergence of homogenization in case needed
call IO_error(600) !******************************************************************************************************
return ! calculate average quantities per ip and post results
endif !$OMP PARALLEL DO
do e=el_start,el_end
if (debugger) then if(CPFEM_in /= 0_pInt) then
!$OMP CRITICAL (write2out) ip_start = CPFEM_in
write (6,*) msg ip_end = CPFEM_in
write (6,*) 'single crystallite convergence reached.',cp_en,CPFEM_in,grain else
write (6,'(a,/,3(3(f12.7,x)/))') 'Lp',CPFEM_Lp(1:3,:,grain,CPFEM_in,cp_en) ip_start = 1
write (6,'(a,/,3(3(f12.7,x)/))') 'Fp (new)',CPFEM_Fp_new(1:3,:,grain,CPFEM_in,cp_en) ip_end = FE_Nips(mesh_element(2,e)) ! iterate over all IPs of this element's type
write (6,'(a,/,3(4(f9.3,x)/))') 'state (new)/ MPa',constitutive_state_new(:,grain,CPFEM_in,cp_en)/1e6_pReal endif
write (6,'(a,/,3(3(f9.3,x)/))') 'P / MPa',PK1/1e6_pReal do i=ip_start,ip_end
write (6,'(a,/,9(9(f9.3,x)/))') 'dP/dF / GPa',dPdF/1e9_pReal do g=1,texture_Ngrains(mesh_element(4,e))
!$OMP END CRITICAL (write2out) volfrac = constitutive_matVolFrac(g,i,e)*constitutive_texVolFrac(g,i,e)
endif CPFEM_PK1_bar(:,:,i,e) = CPFEM_PK1_bar(:,:,i,e) + volfrac * CPFEM_PK1(:,:,g,i,e)
if (updateJaco) CPFEM_dPdF_bar(:,:,:,:,i,e) = &
volfrac = constitutive_matVolFrac(grain,CPFEM_in,cp_en)*constitutive_texVolFrac(grain,CPFEM_in,cp_en) CPFEM_dPdF_bar(:,:,:,:,i,e) + volfrac * CPFEM_dPdF(:,:,:,:,g,i,e) ! add up crystallite stiffnesses
CPFEM_PK1_bar(:,:,CPFEM_in,cp_en) = CPFEM_PK1_bar(:,:,CPFEM_in,cp_en) + volfrac*PK1 ! (may have "holes" corresponding
if (updateJaco) CPFEM_dPdF_bar(:,:,:,:,CPFEM_in,cp_en) = & ! to former avg tangent)
CPFEM_dPdF_bar(:,:,:,:,CPFEM_in,cp_en) + volfrac*dPdF ! add up crystallite stiffnesses ! update results plotted in MENTAT
! (may have "holes" corresponding call math_pDecomposition(CPFEM_Fe1(:,:,g,i,e),U,R,error) ! polar decomposition
! to former avg tangent) if (error) then
! !$OMP CRITICAL (write2out)
! update results plotted in MENTAT write(6,*) 'polar decomposition of', CPFEM_Fe1(:,:,g,i,e)
call math_pDecomposition(Fe1,U,R,error) ! polar decomposition write(6,*) 'Grain: ',g
if (error) then write(6,*) 'Integration point: ',i
!$OMP CRITICAL (write2out) write(6,*) 'Element: ',mesh_element(1,e)
write(6,*) 'polar decomposition of', Fe1 !$OMP END CRITICAL (write2out)
write(6,*) 'Grain: ',grain call IO_error(650)
write(6,*) 'Integration point: ',CPFEM_in return
write(6,*) 'Element: ',mesh_element(1,cp_en) endif
!$OMP END CRITICAL (write2out) CPFEM_results(1:3,g,i,e) = math_RtoEuler(transpose(R))*inDeg ! orientation
call IO_error(650) CPFEM_results(4 ,g,i,e) = volfrac ! volume fraction of orientation
return end do
endif end do
CPFEM_results(1:3,grain,CPFEM_in,cp_en) = math_RtoEuler(transpose(R))*inDeg ! orientation end do
CPFEM_results(4 ,grain,CPFEM_in,cp_en) = volfrac ! volume fraction of orientation !$OMP END PARALLEL DO
enddo ! grain !
! return
return !
! END SUBROUTINE
END SUBROUTINE !
! !
END MODULE !********************************************************************
!############################################################## ! Initialize crystallite
!********************************************************************
subroutine crystallite_init()
use mesh, only: mesh_maxNips,mesh_NcpElems
use constitutive, only: constitutive_maxNgrains
implicit none
allocate(crystallite_converged(constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)); crystallite_converged = .false.
!
! *** Output to MARC output file ***
!$OMP CRITICAL (write2out)
write(6,*)
write(6,*) 'crystallite Initialization'
write(6,*)
write(6,*) 'crystallite_converged: ', shape(crystallite_converged)
write(6,*)
call flush(6)
!$OMP END CRITICAL (write2out)
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
el_start,& ! first element in element loop
el_end,& ! last element in element loop
CPFEM_in) ! IP number
!
use prec, only: pReal,pInt,pert_Fg,subStepMin, nCutback
use debug
use constitutive
use mesh, only: mesh_element, FE_Nips
use math
use IO, only: IO_error
! use CPFEM
implicit none
!
logical updateJaco, JacoOK
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_maxNstatevars) :: state_pert
integer(pInt) el_start, el_end, CPFEM_in, ip_start, ip_end, g, i, e, k, l, iOuter
!
crystallite_converged=.true.
!$OMP PARALLEL DO
do e=el_start,el_end
if(CPFEM_in /= 0_pInt) then
ip_start = CPFEM_in
ip_end = CPFEM_in
else
ip_start = 1
ip_end = FE_Nips(mesh_element(2,e)) ! iterate over all IPs of this element's type
endif
do i=ip_start,ip_end
do g=1,texture_Ngrains(mesh_element(4,e))
crystallite_converged(g,i,e)=.false.
end do
end do
end do
!$OMP END PARALLEL DO
constitutive_state_new=constitutive_state_old
CPFEM_Lp_new = CPFEM_Lp_old
iOuter = 0_pInt
do while(any(crystallite_converged(:,:,el_start:el_end))==.false.)
!$OMP PARALLEL DO
do e=el_start,el_end
if(CPFEM_in /= 0_pInt) then
ip_start = CPFEM_in
ip_end = CPFEM_in
else
ip_start = 1
ip_end = FE_Nips(mesh_element(2,e)) ! iterate over all IPs of this element's type
endif
do i=ip_start,ip_end
do g=1,texture_Ngrains(mesh_element(4,e))
if(.not.crystallite_converged(g,i,e))&
call IntegrateStress(CPFEM_Tstar_v(:,g,i,e), CPFEM_PK1(:,:,g,i,e), CPFEM_ffn1(:,:,g,i,e),&
CPFEM_Fp_new(:,:,g,i,e), CPFEM_Fe1(:,:,g,i,e), CPFEM_Lp_new(:,:,g,i,e),&
constitutive_state_new(:,g,i,e), dt, g, i, e)
end do
end do
end do
!$OMP END PARALLEL DO
!$OMP PARALLEL DO
do e=el_start,el_end
if(CPFEM_in /= 0_pInt) then
ip_start = CPFEM_in
ip_end = CPFEM_in
else
ip_start = 1
ip_end = FE_Nips(mesh_element(2,e)) ! iterate over all IPs of this element's type
endif
do i=ip_start,ip_end
do g=1,texture_Ngrains(mesh_element(4,e))
if(.not.crystallite_converged(g,i,e))&
call UpdateState(CPFEM_Tstar_v(:,g,i,e),constitutive_state_new(:,g,i,e),dt,g,i,e)
end do
end do
end do
!$OMP END PARALLEL DO
iOuter = iOuter + 1_pInt
if (iOuter==Nouter) then
!$OMP CRITICAL (write2out)
write (6,*) 'Terminated outer loop at el,ip,grain',e,i,g
!$OMP CRITICAL (out)
debug_OuterLoopDistribution(iOuter) = debug_OuterLoopDistribution(iOuter)+1
!$OMP END CRITICAL (out)
call IO_error(600)
!$OMP END CRITICAL (write2out)
endif
end do
!$OMP CRITICAL (out)
debug_OuterLoopDistribution(iOuter) = debug_OuterLoopDistribution(iOuter)+1
!$OMP END CRITICAL (out)
if (wantsConstitutiveResults) then ! get the post_results upon request
!$OMP PARALLEL DO
do e=el_start,el_end
if(CPFEM_in /= 0_pInt) then
ip_start = CPFEM_in
ip_end = CPFEM_in
else
ip_start = 1
ip_end = FE_Nips(mesh_element(2,e)) ! iterate over all IPs of this element's type
endif
do i=ip_start,ip_end
do g=1,texture_Ngrains(mesh_element(4,e))
CPFEM_results(CPFEM_Nresults+1:CPFEM_Nresults+constitutive_Nresults(g,i,e),g,i,e) =&
constitutive_post_results(CPFEM_Tstar_v(:,g,i,e),constitutive_state_new(:,g,i,e),&
CPFEM_Temperature(i,e),dt,g,i,e)
end do
end do
end do
!$OMP END PARALLEL DO
endif
!
!***** Calculate Jacobian *****
if(updateJaco) then
if (debugger) then
!$OMP CRITICAL (write2out)
write (6,*) 'Jacobian calc'
!$OMP END CRITICAL (write2out)
endif
! crystallite_converged=.false.
!$OMP PARALLEL DO
do e=el_start,el_end
if(CPFEM_in /= 0_pInt) then
ip_start = CPFEM_in
ip_end = CPFEM_in
else
ip_start = 1
ip_end = FE_Nips(mesh_element(2,e)) ! iterate over all IPs of this element's type
endif
do i=ip_start,ip_end
do g=1,texture_Ngrains(mesh_element(4,e))
do k=1,3
do l=1,3
crystallite_converged(g,i,e)=.false.
JacoOK=.true.
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
state_pert = constitutive_state_new(:,g,i,e) ! initial guess from end of time step
iOuter=0_pInt
do while(.not.crystallite_converged(g,i,e))
call IntegrateStress(Tstar_v, P_pert, Fg_pert, Fp_pert, Fe_pert, Lp_pert, state_pert, dt, g, i, e)
call UpdateState(Tstar_v,state_pert,dt,g,i,e)
iOuter = iOuter + 1_pInt
if (iOuter==Nouter) then
JacoOK=.false.
exit
endif
end do
!$OMP CRITICAL (out)
debug_OuterLoopDistribution(iOuter) = debug_OuterLoopDistribution(iOuter)+1
!$OMP END CRITICAL (out)
if (JacoOK) &
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
! otherwise leave component unchanged
end do
end do
end do
end do
end do
!$OMP END PARALLEL DO
endif
!
return
!
end subroutine
!
!********************************************************************
! Update the state for a single component
!********************************************************************
subroutine UpdateState(&
Tstar_v,& ! stress
state,& ! state
dt,& ! time increment
g,& ! grain number
i,& ! integration point number
e& ! element number
)
use prec, only: pReal,pInt,reltol_Outer
use constitutive, only: constitutive_dotState, constitutive_state_old, constitutive_Nstatevars
! use CPFEM, only: CPFEM_Temperature
!
integer(pInt) g, i, e
real(pReal), dimension(6) :: Tstar_v
real(pReal), dimension(constitutive_Nstatevars(g, i, e)) :: state, ROuter
real(pReal) dt
!
ROuter = state - constitutive_state_old(:,g,i,e) - &
dt*constitutive_dotState(Tstar_v,state,CPFEM_Temperature(i,e),&
g,i,e) ! residuum from evolution of microstructure
state = state - ROuter ! update of microstructure
if (maxval(abs(ROuter/state),state /= 0.0_pReal) < reltol_Outer) crystallite_converged(g,i,e) = .true.
!
return
!
end subroutine
!
!
!********************************************************************
! Calculates the stress for a single component
!********************************************************************
!***********************************************************************
!*** calculation of stress (P), stiffness (dPdF), ***
!*** and announcment of any ***
!*** acceleration of the Newton-Raphson correction ***
!***********************************************************************
subroutine IntegrateStress(&
Tstar_v,& ! Stress vector
P,& ! first PK stress
Fg_new,& ! new global deformation gradient
Fp_new,& ! new plastic deformation gradient
Fe_new,& ! new "elastic" deformation gradient
Lp,& ! plastic velocity gradient
state_new,& ! new state variable array
dt,& ! time increment
g,& ! grain number
i,& ! integration point number
e) ! element number
! post_results,& ! plot results from constitutive model
! Fp_new,& ! new plastic deformation gradient
! updateJaco,& ! update of Jacobian required
! Temperature,& ! temperature of crystallite
! Fg_old,& ! old global deformation gradient
! Fp_old,& ! old plastic deformation gradient
! state_old) ! old state variable array
!
use prec, only: pReal,pInt,pert_Fg,subStepMin, nCutback
use debug
use constitutive, only: constitutive_Nstatevars,constitutive_Nresults,constitutive_state_old
use math
! use CPFEM
!
implicit none
!
character(len=128) 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
real(pReal), dimension(constitutive_Nstatevars(g,i,e)) :: state_new
! real(pReal), dimension(constitutive_Nstatevars(g,i,e)) :: state_current
!
! debugger= e==1.and.i==1
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)
! state_current = state_new
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
! state_current = state_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 current / MPa',state_current/1e6_pReal
write (6,'(a,/,3(4(f9.3,x)/))') 'state new / MPa',state_new/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,state_new,dt_aim,e,i,g,Temperature,Fg_aim,Fp_current)
!
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
! if (debugger) then
!$OMP CRITICAL (write2out)
write (6,*) '>>>>>>>>>>>>>>>>>>>> cutback <<<<<<<<<<<<<<<<<<<<<<'
write (6,*) 'Element, Ip:', e, i
write (6,*) msg
!$OMP END CRITICAL (write2out)
! endif
!
endif
enddo ! potential substepping
!
!$OMP CRITICAL (cutback)
debug_cutbackDistribution(min(nCutback,maxCutbacks)+1) = debug_cutbackDistribution(min(nCutback,maxCutbacks)+1)+1
!$OMP END CRITICAL (cutback)
!
! debugger = .false.
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,& ! 1nd PK stress (taken as initial guess if /= 0)
state,& ! current microstructure at end of time inc (taken as guess if /= 0)
dt,& ! time increment
cp_en,& ! element number
ip,& ! integration point number
grain,& ! grain number
Temperature,& ! temperature
Fg_new,& ! new total def gradient
Fp_old) ! former plastic def gradient
! state_current) ! former microstructure
use prec
use debug
use mesh, only: mesh_element
use constitutive, only: constitutive_Nstatevars,constitutive_Microstructure,&
constitutive_homogenizedC,constitutive_LpAndItsTangent
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
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 !,Tstar
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
real(pReal), dimension(constitutive_Nstatevars(grain, ip, cp_en)) :: state
!
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)))
!
! if (all(state == 0.0_pReal)) state = state_current ! former state guessed, if none specified
! iOuter = 0_pInt ! outer counter
!
if (debugger) then
!$OMP CRITICAL (write2out)
write (6,'(a,/,3(3(f12.7,x)/))') 'Fg to be calculated',Fg_new
!$OMP END CRITICAL (write2out)
endif
!
!Outer: do ! outer iteration: State
! iOuter = iOuter+1
! if (debugger) then
!!$OMP CRITICAL (write2out)
! write (6,'(a,i3)') '---outer ',iOuter
! write (6,'(a,/,3(4(f9.3,x)/))') 'state old / MPa',state_old/1e6_pReal
! write (6,'(a,/,3(4(f9.3,x)/))') 'state / MPa',state/1e6_pReal
! write (6,'(a,/,3(3(f12.7,x)/))') 'Lpguess',Lpguess(1:3,:)
!!$OMP END CRITICAL (write2out)
! endif
!
! if (iOuter > nOuter) then
! msg = 'limit Outer iteration'
!!$OMP CRITICAL (out)
! debug_OuterLoopDistribution(nOuter) = debug_OuterLoopDistribution(nOuter)+1
!!$OMP END CRITICAL (out)
! return
! endif
call constitutive_Microstructure(state,Temperature,grain,ip,cp_en)
C_66 = constitutive_HomogenizedC(state, 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 (debugger) then
!!$OMP CRITICAL (write2out)
! write (6,'(a,i3)') 'inner ',iInner
! if (iInner < 3) then
! write (6,'(a,/,3(3(f12.7,x)/))') 'Lpguess',Lpguess(1:3,:)
! endif
!!$OMP END CRITICAL (write2out)
! endif
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'
!$OMP CRITICAL (in)
debug_InnerLoopDistribution(nInner) = debug_InnerLoopDistribution(nInner)+1
!$OMP END CRITICAL (in)
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))
! Tstar = math_Mandel6to33(Tstar_v)
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 constitutive_LpAndItsTangent(Lp,dLp, &
Tstar_v,state,Temperature,grain,ip,cp_en)
!
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 / 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 (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)
! ROuter = state - state_old - &
! dt*constitutive_dotState(Tstar_v,state,Temperature,&
! grain,ip,cp_en) ! residuum from evolution of microstructure
! state = state - ROuter ! update of microstructure
!
! if (iOuter==nOuter) then
!!$OMP CRITICAL (write2out)
! write (6,*) 'Terminated outer loop at el,ip,grain',cp_en,ip,grain
!!$OMP END CRITICAL (write2out)
! exit Outer
! endif
! if (maxval(abs(Router/state),state /= 0.0_pReal) < reltol_Outer) exit Outer
! enddo Outer
!
!!$OMP CRITICAL (out)
! debug_OuterLoopDistribution(iOuter) = debug_OuterLoopDistribution(iOuter)+1
!!$OMP END CRITICAL (out)
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
!
! if (wantsConstitutiveResults) then ! get the post_results upon request
! results = 0.0_pReal
! results = constitutive_post_results(Tstar_v,state,Temperature,dt,grain,ip,cp_en)
! 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(Tstar,transpose(invFp_new))) ! first PK stress
P = math_mul33x33(Fe_new,math_mul33x33(math_Mandel6to33(Tstar_v),transpose(invFp_new))) ! first PK stress
return
!
END SUBROUTINE
!
END MODULE
!##############################################################

398
trunk/crystallite.f90
View File

@ -1,398 +0,0 @@
!##############################################################
MODULE crystallite
!##############################################################
! *** Solution at single crystallite level ***
!
CONTAINS
!
!
!********************************************************************
! Calculates the stress for a single component
!********************************************************************
!***********************************************************************
!*** calculation of stress (P), stiffness (dPdF), ***
!*** and announcment of any ***
!*** acceleration of the Newton-Raphson correction ***
!***********************************************************************
subroutine SingleCrystallite(&
msg,& ! return message
P,& ! first PK stress
dPdF,& ! consistent tangent
post_results,& ! plot results from constitutive model
Lp,& ! plastic velocity gradient
Fp_new,& ! new plastic deformation gradient
Fe_new,& ! new "elastic" deformation gradient
state_new,& ! new state variable array
dt,& ! time increment
cp_en,& ! element number
ip,& ! integration point number
grain,& ! grain number
updateJaco,& ! update of Jacobian required
Temperature,& ! temperature of crystallite
Fg_new,& ! new global deformation gradient
Fg_old,& ! old global deformation gradient
Fp_old,& ! old plastic deformation gradient
state_old) ! old state variable array
!
use prec, only: pReal,pInt,pert_Fg,subStepMin, nCutback
use debug
use constitutive, only: constitutive_Nstatevars,constitutive_Nresults
use mesh, only: mesh_element
use math
implicit none
!
character(len=*) msg
logical updateJaco,error,success,guessNew
integer(pInt) cp_en,ip,grain,k,l, nCutbacks, maxCutbacks
real(pReal) Temperature
real(pReal) dt,dt_aim,subFrac,subStep,det
real(pReal), dimension(3,3) :: Lp,Lp_interpolated,Lp_pert,inv
real(pReal), dimension(3,3) :: Fg_old,Fg_current,Fg_new,Fg_pert,Fg_aim,deltaFg
real(pReal), dimension(3,3) :: Fp_old,Fp_current,Fp_new,Fp_pert
real(pReal), dimension(3,3) :: Fe_current,Fe_new,Fe_pert
real(pReal), dimension(3,3) :: P,P_pert
real(pReal), dimension(3,3,3,3) :: dPdF
real(pReal), dimension(constitutive_Nstatevars(grain,ip,cp_en)) :: state_old,state_new
real(pReal), dimension(constitutive_Nstatevars(grain,ip,cp_en)) :: state_current,state_bestguess,state_pert
real(pReal), dimension(constitutive_Nresults(grain,ip,cp_en)) :: post_results
deltaFg = Fg_new - Fg_old
subFrac = 0.0_pReal
subStep = 1.0_pReal
nCutbacks = 0_pInt
maxCutbacks = 0_pInt
Fg_current = Fg_old ! initialize to start of inc
Fp_current = Fp_old
call math_invert3x3(Fp_old,inv,det,error)
Fe_current = math_mul33x33(Fg_old,inv)
state_current = state_old
success = .false. ! pretend cutback
dt_aim = 0.0_pReal ! prevent initial Lp interpolation
! 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
state_current = state_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 current / MPa',state_current/1e6_pReal
write (6,'(a,/,3(4(f9.3,x)/))') 'state new / MPa',state_new/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,P,state_new,post_results, & ! def gradients and PK2 at end of time step
maxval(abs(Fg_aim-Fg_new)) < relevantStrain, & ! post results only if asking for final values
dt_aim,cp_en,ip,grain,Temperature,Fg_aim,Fp_current,state_current)
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
if (debugger) then
!$OMP CRITICAL (write2out)
write (6,*) '>>>>>>>>>>>>>>>>>>>> cutback <<<<<<<<<<<<<<<<<<<<<<'
!$OMP END CRITICAL (write2out)
endif
endif
enddo ! potential substepping
!
!$OMP CRITICAL (cutback)
debug_cutbackDistribution(min(nCutback,maxCutbacks)+1) = debug_cutbackDistribution(min(nCutback,maxCutbacks)+1)+1
!$OMP END CRITICAL (cutback)
!
if (msg /= 'ok') return ! solution not reached --> report back
if (updateJaco) then ! consistent tangent using
if (debugger) then
!$OMP CRITICAL (write2out)
write (6,*) 'Jacobian calc'
!$OMP END CRITICAL (write2out)
endif
do k=1,3
do l=1,3
Fg_pert = Fg_new ! initialize perturbed Fg
Fg_pert(k,l) = Fg_pert(k,l) + pert_Fg ! perturb single component
Lp_pert = Lp
state_pert = state_new ! initial guess from end of time step
call TimeIntegration(msg,Lp_pert,Fp_pert,Fe_pert,P_pert,state_pert,post_results,.false., & ! def gradients and PK2 at end of time step
dt_aim,cp_en,ip,grain,Temperature,Fg_pert,Fp_current,state_current)
if (msg == 'ok') &
dPdF(:,:,k,l) = (P_pert-P)/pert_Fg ! constructing tangent dP_ij/dFg_kl only if valid forward difference
! otherwise leave component unchanged
enddo
enddo
endif
!
msg = 'ok' ! a new consistent tangent was computed even if msg was not ok for all components
!
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
P,& ! 1nd PK stress (taken as initial guess if /= 0)
state,& ! current microstructure at end of time inc (taken as guess if /= 0)
results,& ! post results from constitutive
wantsConstitutiveResults,& ! its flag
!
dt,& ! time increment
cp_en,& ! element number
ip,& ! integration point number
grain,& ! grain number
Temperature,& ! temperature
Fg_new,& ! new total def gradient
Fp_old,& ! former plastic def gradient
state_old) ! former microstructure
use prec
use debug
use mesh, only: mesh_element
use constitutive, only: constitutive_Nstatevars,&
constitutive_homogenizedC,constitutive_dotState,constitutive_LpAndItsTangent,&
constitutive_Nresults,constitutive_Microstructure,constitutive_post_results
use math
use IO
implicit none
!
character(len=*) msg
logical failed,wantsConstitutiveResults
integer(pInt) cp_en, ip, grain
integer(pInt) iOuter,iInner,dummy, i,j,k,l,m,n
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,Tstar
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
real(pReal), dimension(constitutive_Nstatevars(grain, ip, cp_en)) :: state_old,state,ROuter
real(pReal), dimension(constitutive_Nresults(grain,ip,cp_en)) :: results
!
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)))
!
if (all(state == 0.0_pReal)) state = state_old ! former state guessed, if none specified
iOuter = 0_pInt ! outer counter
!
if (debugger) then
!$OMP CRITICAL (write2out)
write (6,'(a,/,3(3(f12.7,x)/))') 'Fg to be calculated',Fg_new
!$OMP END CRITICAL (write2out)
endif
!
Outer: do ! outer iteration: State
iOuter = iOuter+1
if (debugger) then
!$OMP CRITICAL (write2out)
write (6,'(a,i3)') '---outer ',iOuter
write (6,'(a,/,3(4(f9.3,x)/))') 'state old / MPa',state_old/1e6_pReal
write (6,'(a,/,3(4(f9.3,x)/))') 'state / MPa',state/1e6_pReal
write (6,'(a,/,3(3(f12.7,x)/))') 'Lpguess',Lpguess(1:3,:)
!$OMP END CRITICAL (write2out)
endif
if (iOuter > nOuter) then
msg = 'limit Outer iteration'
!$OMP CRITICAL (out)
debug_OuterLoopDistribution(nOuter) = debug_OuterLoopDistribution(nOuter)+1
!$OMP END CRITICAL (out)
return
endif
call constitutive_Microstructure(state,Temperature,grain,ip,cp_en)
C_66 = constitutive_HomogenizedC(state, 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 (debugger) then
!$OMP CRITICAL (write2out)
write (6,'(a,i3)') 'inner ',iInner
if (wantsConstitutiveResults .and. iOuter == 1 .and. iInner < 3) then
write (6,'(a,/,3(3(f12.7,x)/))') 'Lpguess',Lpguess(1:3,:)
endif
!$OMP END CRITICAL (write2out)
endif
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'
!$OMP CRITICAL (in)
debug_InnerLoopDistribution(nInner) = debug_InnerLoopDistribution(nInner)+1
!$OMP END CRITICAL (in)
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))
Tstar = math_Mandel6to33(Tstar_v)
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 constitutive_LpAndItsTangent(Lp,dLp, &
Tstar_v,state,Temperature,grain,ip,cp_en)
!
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 / 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 (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)
ROuter = state - state_old - &
dt*constitutive_dotState(Tstar_v,state,Temperature,&
grain,ip,cp_en) ! residuum from evolution of microstructure
state = state - ROuter ! update of microstructure
if (iOuter==nOuter) then
!!$OMP CRITICAL (write2out)
write (6,*) 'Terminated outer loop at el,ip,grain',cp_en,ip,grain
!!$OMP END CRITICAL (write2out)
exit Outer
endif
if (maxval(abs(Router/state),state /= 0.0_pReal) < reltol_Outer) exit Outer
enddo Outer
!
!$OMP CRITICAL (out)
debug_OuterLoopDistribution(iOuter) = debug_OuterLoopDistribution(iOuter)+1
!$OMP END CRITICAL (out)
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
!
if (wantsConstitutiveResults) then ! get the post_results upon request
results = 0.0_pReal
results = constitutive_post_results(Tstar_v,state,Temperature,dt,grain,ip,cp_en)
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(Tstar,transpose(invFp_new))) ! first PK stress
return
!
END SUBROUTINE
!
!
END MODULE
!##############################################################

332
trunk/mpie_cpfem_marc2005r3.f90
View File

@ -1,332 +0,0 @@
!********************************************************************
! Material subroutine for MSC.Marc Version 0.1
!
! written by F. Roters, P. Eisenlohr, L. Hantcherli, W.A. Counts
! MPI fuer Eisenforschung, Duesseldorf
!
! last modified: 08.11.2007
!********************************************************************
! Usage:
! - choose material as hypela2
! - set statevariable 2 to index of material
! - set statevariable 3 to index of texture
! - choose output of user variables if desired
! - make sure the file "mattex.mpie" exists in the working
! directory
! - use nonsymmetric option for solver (e.g. direct
! profile or multifrontal sparse, the latter seems
! to be faster!)
!********************************************************************
! Marc subroutines used:
! - hypela2
! - plotv
! - quit
!********************************************************************
! Marc common blocks included:
! - concom: lovl, ncycle, inc, incsub
! - creeps: timinc
!********************************************************************
!
include "prec.f90" ! uses nothing else
include "debug.f90" ! uses prec
include "math.f90" ! uses prec
include "IO.f90" ! uses prec, debug, math
include "FEsolving.f90" ! uses prec, IO
include "mesh.f90" ! uses prec, IO, math, FEsolving
include "lattice.f90" ! uses prec, math
include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug
include "crystallite.f90" ! uses prec, debug, constitutive, mesh, math, IO
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
!********************************************************************
!
! ************* user subroutine for defining material behavior **************
!
!
! CAUTION : Due to calculation of the Deformation gradients, Stretch Tensors and
! Rotation tensors at previous and current states, the analysis can be
! computationally expensive. Please use the user subroutine -> hypela
! if these kinematic quantities are not needed in the constitutive model
!
!
! IMPORTANT NOTES :
!
! (1) F,R,U are only available for continuum and membrane elements (not for
! shells and beams).
!
! (2) For total Lagrangian formulation use the -> 'Elasticity,1' card(=
! total Lagrange with large disp) in the parameter section of input deck.
! For updated Lagrangian formulation use the -> 'Plasticity,3' card(=
! update+finite+large disp+constant d) in the parameter section of
! 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 ! Marc common blocks are in fixed format so they have to be pasted in here
!
! Marc common blocks are in fixed format so they have to be pasted in here
! Beware of changes in newer Marc versions -- these are from 2005r3
! concom is needed for inc, subinc, ncycle, lovl
! include 'concom'
common/concom/ &
iacous, iasmbl, iautth, ibear, icompl, iconj, icreep, ideva(50), idyn, idynt,&
ielas, ielcma, ielect, iform, ifour, iharm, ihcps, iheat, iheatt, ihresp,&
ijoule, ilem, ilnmom, iloren, inc, incext, incsub, ipass, iplres, ipois,&
ipoist, irpflo, ismall, ismalt, isoil, ispect, ispnow, istore, iswep, ithcrp,&
itherm, iupblg, iupdat, jacflg, jel, jparks, largst, lfond, loadup, loaduq,&
lodcor, lovl, lsub, magnet, ncycle, newtnt, newton, noshr, linear, ivscpl,&
icrpim, iradrt, ipshft, itshr, iangin, iupmdr, iconjf, jincfl, jpermg, jhour,&
isolvr, jritz, jtable, jshell, jdoubl, jform, jcentr, imini, kautth, iautof,&
ibukty, iassum, icnstd, icnstt, kmakmas, imethvp,iradrte,iradrtp, iupdate,iupdatp,&
ncycnt, marmen ,idynme, ihavca, ispf, kmini, imixed, largtt, kdoela, iautofg,&
ipshftp,idntrc, ipore, jtablm, jtablc, isnecma,itrnspo,imsdif, jtrnspo,mcnear,&
imech, imecht, ielcmat, ielectt,magnett, imsdift,noplas, jtabls, jactch, jtablth,&
kgmsto ,jpzo, ifricsh, iremkin,iremfor, ishearp,jspf, machining, jlshell,icompsol,&
iupblgfo,jcondir,nstcrp, nactive,ipassref, nstspnt,ibeart,icheckmpc, noline, icuring,&
ishrink,ioffsflg,isetoff, iharmt, inc_incdat, iautspc,ibrake
! creeps is needed for timinc (time increment)
! include 'creeps'
common/marc_creeps/ &
cptim,timinc,timinc_p,timinc_s,timincm,timinc_a,timinc_b,creept(33),icptim,icfte,icfst,&
icfeq,icftm,icetem,mcreep,jcreep,icpa,icftmp,icfstr,icfqcp,icfcpm,icrppr,icrcha,icpb,iicpmt,iicpa
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 :
!
! call scla(un1,0.d0,itel,itel,1)
! do 3 k=1,3
! do 2 i=1,3
! do 1 j=1,3
! un1(i,j)=un1(i,j)+dsqrt(strechn1(k))*eigvn1(i,k)*eigvn1(j,k)
!1 continue
!2 continue
!3 continue
!
use prec, only: pReal,pInt, ijaco
use FEsolving
use CPFEM, only: CPFEM_general
use math, only: invnrmMandel
implicit real(pReal) (a-h,o-z)
integer(pInt) computationMode
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)
! Marc common blocks are in fixed format so they have to be pasted in here
! Beware of changes in newer Marc versions -- these are from 2005r3
! concom is needed for inc, subinc, ncycle, lovl
! include 'concom'
common/marc_concom/ &
iacous, iasmbl, iautth, ibear, icompl, iconj, icreep, ideva(50), idyn, idynt,&
ielas, ielcma, ielect, iform, ifour, iharm, ihcps, iheat, iheatt, ihresp,&
ijoule, ilem, ilnmom, iloren, inc, incext, incsub, ipass, iplres, ipois,&
ipoist, irpflo, ismall, ismalt, isoil, ispect, ispnow, istore, iswep, ithcrp,&
itherm, iupblg, iupdat, jacflg, jel, jparks, largst, lfond, loadup, loaduq,&
lodcor, lovl, lsub, magnet, ncycle, newtnt, newton, noshr, linear, ivscpl,&
icrpim, iradrt, ipshft, itshr, iangin, iupmdr, iconjf, jincfl, jpermg, jhour,&
isolvr, jritz, jtable, jshell, jdoubl, jform, jcentr, imini, kautth, iautof,&
ibukty, iassum, icnstd, icnstt, kmakmas, imethvp,iradrte,iradrtp, iupdate,iupdatp,&
ncycnt, marmen ,idynme, ihavca, ispf, kmini, imixed, largtt, kdoela, iautofg,&
ipshftp,idntrc, ipore, jtablm, jtablc, isnecma,itrnspo,imsdif, jtrnspo,mcnear,&
imech, imecht, ielcmat, ielectt,magnett, imsdift,noplas, jtabls, jactch, jtablth,&
kgmsto ,jpzo, ifricsh, iremkin,iremfor, ishearp,jspf, machining, jlshell,icompsol,&
iupblgfo,jcondir,nstcrp, nactive,ipassref, nstspnt,ibeart,icheckmpc, noline, icuring,&
ishrink,ioffsflg,isetoff, ioffsetm,iharmt, inc_incdat,iautspc,ibrake, icbush ,istream_input,&
iprsinp,ivlsinp,ifirst_time,ipin_m,jgnstr_glb, imarc_return,iqvcinp,nqvceid,istpnx,imicro1
! creeps is needed for timinc (time increment)
! include 'creeps'
common/marc_creeps/ &
cptim,timinc,timinc_p,timinc_s,timincm,timinc_a,timinc_b,creept(33),icptim,icfte,icfst,&
icfeq,icftm,icetem,mcreep,jcreep,icpa,icftmp,icfstr,icfqcp,icfcpm,icrppr,icrcha,icpb,iicpmt,iicpa
if (inc == 0) then
cycleCounter = 0
else
if (theCycle > ncycle) cycleCounter = 0 ! reset counter for each cutback
if (theCycle /= ncycle .or. theLovl /= lovl) cycleCounter = cycleCounter+1 ! ping pong
endif
if (cptim > theTime .or. theInc /= inc) then ! reached convergence
lastIncConverged = .true.
outdatedByNewInc = .true.
endif
if (mod(cycleCounter,2) /= 0) computationMode = 4 ! recycle
if (mod(cycleCounter,4) == 2) computationMode = 3 ! collect
if (mod(cycleCounter,4) == 0) computationMode = 2 ! compute
if (computationMode == 4 .and. ncycle == 0 .and. .not. lastIncConverged) &
computationMode = 6 ! recycle but restore known good consistent tangent
if (computationMode == 4 .and. lastIncConverged) then
computationMode = 5 ! recycle and record former consistent tangent
lastIncConverged = .false.
endif
if (computationMode == 2 .and. outdatedByNewInc) then
computationMode = 1 ! compute and age former results
outdatedByNewInc = .false.
endif
theTime = cptim ! record current starting time
theInc = inc ! record current increment number
theCycle = ncycle ! record current cycle count
theLovl = lovl ! record current lovl
call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,n(1),nn,s,mod(theCycle,2_pInt*ijaco)==0,d,ngens)
! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13
! Marc: 11, 22, 33, 12, 23, 13
forall(i=1:ngens) d(1:ngens,i) = invnrmMandel(i)*d(1:ngens,i)*invnrmMandel(1:ngens)
s(1:ngens) = s(1:ngens)*invnrmMandel(1:ngens)
if(symmetricSolver) d(1:ngens,1:ngens) = 0.5_pReal*(d(1:ngens,1:ngens)+transpose(d(1:ngens,1:ngens)))
return
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
!********************************************************************
!
! 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
!
!********************************************************************
use prec, only: pReal,pInt
use CPFEM, only: CPFEM_results, CPFEM_Nresults
use constitutive, only: constitutive_maxNresults
use mesh, only: mesh_FEasCP
implicit none
!
real(pReal) s(*),etot(*),eplas(*),ecreep(*),sp(*)
real(pReal) v, t(*)
integer(pInt) m, nn, layer, ndi, nshear, jpltcd
!
! assign result variable
v=CPFEM_results(mod(jpltcd-1_pInt, CPFEM_Nresults+constitutive_maxNresults)+1_pInt,&
(jpltcd-1_pInt)/(CPFEM_Nresults+constitutive_maxNresults)+1_pInt,&
nn, mesh_FEasCP('elem', m))
return
END SUBROUTINE
!
!
! subroutine utimestep(timestep,timestepold,icall,time,timeloadcase)
!********************************************************************
! This routine modifies the addaptive time step of Marc
!********************************************************************
! use prec, only: pReal,pInt
! use CPFEM, only : CPFEM_timefactor_max
! implicit none
!
! real(pReal) timestep, timestepold, time,timeloadcase
! integer(pInt) icall
!
! user subroutine for modifying the time step in auto step
!
! timestep : the current time step as suggested by marc
! to be modified in this routine
! timestepold : the current time step before it was modified by marc
! icall : =1 for setting the initial time step
! =2 if this routine is called during an increment
! =3 if this routine is called at the beginning
! of the increment
! time : time at the start of the current increment
! timeloadcase: time period of the current load case
!
! it is in general not recommended to increase the time step
! during the increment.
! this routine is called right after the time step has (possibly)
! been updated by marc.
!
! user coding
! reduce timestep during increment in case mpie_timefactor is too large
! if(icall==2_pInt) then
! if(mpie_timefactor_max>1.25_pReal) then
! timestep=min(timestep,timestepold*0.8_pReal)
! end if
! return
! modify timestep at beginning of new increment
! else if(icall==3_pInt) then
! if(mpie_timefactor_max<=0.8_pReal) then
! timestep=min(timestep,timestepold*1.25_pReal)
! else if (mpie_timefactor_max<=1.0_pReal) then
! timestep=min(timestep,timestepold/mpie_timefactor_max)
! else if (mpie_timefactor_max<=1.25_pReal) then
! timestep=min(timestep,timestepold*1.01_pReal)
! else
! timestep=min(timestep,timestepold*0.8_pReal)
! end if
! end if
! return
! end

4
trunk/mpie_cpfem_marc2007r1.f90
View File

@ -4,7 +4,7 @@
! written by F. Roters, P. Eisenlohr, L. Hantcherli, W.A. Counts ! written by F. Roters, P. Eisenlohr, L. Hantcherli, W.A. Counts
! MPI fuer Eisenforschung, Duesseldorf ! MPI fuer Eisenforschung, Duesseldorf
! !
! last modified: 09.07.2008 ! last modified: 27.11.2008
!******************************************************************** !********************************************************************
! Usage: ! Usage:
! - choose material as hypela2 ! - choose material as hypela2
@ -35,7 +35,7 @@
include "mesh.f90" ! uses prec, IO, math, FEsolving include "mesh.f90" ! uses prec, IO, math, FEsolving
include "lattice.f90" ! uses prec, math include "lattice.f90" ! uses prec, math
include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug
include "crystallite.f90" ! uses prec, debug, constitutive, mesh, math, IO ! include "crystallite.f90" ! uses prec, debug, constitutive, mesh, math, IO
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
! !

4
trunk/mpie_cpfem_marc2007r1_sequential.f90
View File

@ -4,7 +4,7 @@
! written by F. Roters, P. Eisenlohr, L. Hantcherli, W.A. Counts ! written by F. Roters, P. Eisenlohr, L. Hantcherli, W.A. Counts
! MPI fuer Eisenforschung, Duesseldorf ! MPI fuer Eisenforschung, Duesseldorf
! !
! last modified: 09.07.2008 ! last modified: 27.11.2008
!******************************************************************** !********************************************************************
! Usage: ! Usage:
! - choose material as hypela2 ! - choose material as hypela2
@ -35,7 +35,7 @@
include "mesh.f90" ! uses prec, IO, math, FEsolving include "mesh.f90" ! uses prec, IO, math, FEsolving
include "lattice.f90" ! uses prec, math include "lattice.f90" ! uses prec, math
include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug
include "crystallite.f90" ! uses prec, debug, constitutive, mesh, math, IO ! include "crystallite.f90" ! uses prec, debug, constitutive, mesh, math, IO
include "CPFEM_sequential.f90" ! uses prec, math, mesh, constitutive, FEsolving, debug, lattice, IO, crystallite include "CPFEM_sequential.f90" ! uses prec, math, mesh, constitutive, FEsolving, debug, lattice, IO, crystallite
! !

4
trunk/mpie_cpfem_marc2008r1.f90
View File

@ -4,7 +4,7 @@
! written by F. Roters, P. Eisenlohr, L. Hantcherli, W.A. Counts ! written by F. Roters, P. Eisenlohr, L. Hantcherli, W.A. Counts
! MPI fuer Eisenforschung, Duesseldorf ! MPI fuer Eisenforschung, Duesseldorf
! !
! last modified: 28.10.2008 ! last modified: 22.11.2008
!******************************************************************** !********************************************************************
! Usage: ! Usage:
! - choose material as hypela2 ! - choose material as hypela2
@ -36,7 +36,7 @@
include "mesh.f90" ! uses prec, IO, math, FEsolving include "mesh.f90" ! uses prec, IO, math, FEsolving
include "lattice.f90" ! uses prec, math include "lattice.f90" ! uses prec, math
include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug
include "crystallite.f90" ! uses prec, debug, constitutive, mesh, math, IO ! include "crystallite.f90" ! uses prec, debug, constitutive, mesh, math, IO
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
! !
! !

5
trunk/mpie_cpfem_marc2008r1_sequential.f90
View File

@ -4,7 +4,7 @@
! written by F. Roters, P. Eisenlohr, L. Hantcherli, W.A. Counts ! written by F. Roters, P. Eisenlohr, L. Hantcherli, W.A. Counts
! MPI fuer Eisenforschung, Duesseldorf ! MPI fuer Eisenforschung, Duesseldorf
! !
! last modified: 09.07.2008 ! last modified: 27.11.2008
!******************************************************************** !********************************************************************
! Usage: ! Usage:
! - choose material as hypela2 ! - choose material as hypela2
@ -16,6 +16,7 @@
! - use nonsymmetric option for solver (e.g. direct ! - use nonsymmetric option for solver (e.g. direct
! profile or multifrontal sparse, the latter seems ! profile or multifrontal sparse, the latter seems
! to be faster!) ! to be faster!)
! - in case of ddm a symmetric solver has to be used
!******************************************************************** !********************************************************************
! Marc subroutines used: ! Marc subroutines used:
! - hypela2 ! - hypela2
@ -35,7 +36,7 @@
include "mesh.f90" ! uses prec, IO, math, FEsolving include "mesh.f90" ! uses prec, IO, math, FEsolving
include "lattice.f90" ! uses prec, math include "lattice.f90" ! uses prec, math
include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug
include "crystallite.f90" ! uses prec, debug, constitutive, mesh, math, IO ! include "crystallite.f90" ! uses prec, debug, constitutive, mesh, math, IO
include "CPFEM_sequential.f90" ! uses prec, math, mesh, constitutive, FEsolving, debug, lattice, IO, crystallite include "CPFEM_sequential.f90" ! uses prec, math, mesh, constitutive, FEsolving, debug, lattice, IO, crystallite
! !

13
trunk/todo → trunk/todo.txt
View File

@ -1,5 +1,8 @@
Things to be implemented into the code Things to be implemented into the code
# parsing of texture data from mattex file to be done by separate "texture.f90", thus freeing constitutive.f90 from this global task # parsing of texture data from mattex file to be done by separate "texture.f90", thus freeing constitutive.f90 from this global task
# constitutive law to be used is indicated within [material]. constitutive.f90 then assumes a wrapper functionality for the four functions "dor_microstructure", "Lpanditstangent", "postresults" and "microstructure". Within those, a switch case checks for the true constitutive law (a string) and passes control on to the respective subroutine-clone of this law. Technically, the compilation relies on "IFDEF" to include the necessary constitutive laws plus the corresponding entries in each switch case. this should go along with new scheme for multi material definition
# change state variable meaning to (i) homogenization, (ii) microstructure # constitutive law to be used is indicated within [material]. constitutive.f90 then assumes a wrapper functionality for the four functions "dot_microstructure", "Lpanditstangent", "postresults" and "microstructure". Within those, a switch case checks for the true constitutive law (a string) and passes control on to the respective subroutine-clone of this law. Technically, the compilation relies on "IFDEF" to include the necessary constitutive laws plus the corresponding entries in each switch case.
# change state variable meaning to (i) homogenization, (ii) microstructure
# adopt CPFEM_GIA8.f90 to new scheme, rename to CPFEM_RGC.f90
# make OpenMP parallelization work again