diff --git a/trunk/CPFEM.f90 b/trunk/CPFEM.f90 index 38317cf1a..4ad3dd916 100644 --- a/trunk/CPFEM.f90 +++ b/trunk/CPFEM.f90 @@ -9,33 +9,10 @@ ! **************************************************************** ! *** General variables for the material behaviour calculation *** ! **************************************************************** - real(pReal), dimension (:,:), allocatable :: CPFEM_Temperature - real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_ffn_bar !average FFN per IP - real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_ffn !individual FFN per grain - 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 + real(pReal), dimension (:,:,:), allocatable :: CPFEM_cs ! Cauchy stress + real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_dcsdE ! Cauchy stress tangent + real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_dcsdE_knownGood ! known good tangent - 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_calc_done = .false. ! remember whether first IP has already calced the results @@ -47,89 +24,35 @@ !*** allocate the arrays defined in module CPFEM *** !*** and initialize them *** !********************************************************* - SUBROUTINE CPFEM_init(Temperature) -! - use prec - use math, only: math_EulertoR, math_I3, math_identity2nd - use FEsolving, only: parallelExecution - use mesh - use material - use constitutive -! + SUBROUTINE CPFEM_init() + + use prec, only: pInt,pReal + use FEsolving, only: parallelExecution,symmetricSolver,FEsolving_execElem,FEsolving_execIP + use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips,FE_Nips + use material, only: homogenization_maxNgrains + use constitutive, only: constitutive_maxSizePostResults + use crystallite, only: crystallite_Nresults + use homogenization, only: homogenization_maxSizePostResults + implicit none -! - real(pReal) Temperature 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_Fp_new(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Fp_new = 0.0_pReal - allocate(CPFEM_Fp_old(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) - 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_cs(6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_cs = 0.0_pReal + allocate(CPFEM_dcsdE(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dcsde = 0.0_pReal + allocate(CPFEM_dcsdE_knownGood(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dcsde_knownGood = 0.0_pReal - 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 *** !$OMP CRITICAL (write2out) write(6,*) - write(6,*) 'CPFEM Initialization' + write(6,*) '<<<+- cpfem init -+>>>' write(6,*) - write(6,*) 'CPFEM_Temperature: ', shape(CPFEM_Temperature) - write(6,*) 'CPFEM_ffn_bar: ', shape(CPFEM_ffn_bar) - write(6,*) 'CPFEM_ffn: ', shape(CPFEM_ffn) - 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,'(a32,x,6(i5,x))') 'CPFEM_cs: ', shape(CPFEM_cs) + write(6,'(a32,x,6(i5,x))') 'CPFEM_dcsde: ', shape(CPFEM_dcsde) + write(6,'(a32,x,6(i5,x))') 'CPFEM_dcsde_knownGood: ', shape(CPFEM_dcsde_knownGood) write(6,*) write(6,*) 'parallelExecution: ', parallelExecution + write(6,*) 'symmetricSolver: ', symmetricSolver call flush(6) !$OMP END CRITICAL (write2out) return @@ -153,23 +76,27 @@ ! CPFEM_jaco jacobian in Mandel notation ! CPFEM_ngens size of stress strain law !*********************************************************************** - SUBROUTINE CPFEM_general(CPFEM_mode, ffn, ffn1, Temperature, CPFEM_dt,& - CPFEM_en, CPFEM_in, CPFEM_stress, CPFEM_updateJaco, CPFEM_jaco, CPFEM_ngens) +subroutine CPFEM_general(CPFEM_mode, ffn, ffn1, Temperature, CPFEM_dt,& + 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 ! - use prec, only: pReal,pInt + use prec, only: pReal,pInt use FEsolving use debug use math - use mesh, only: mesh_init,mesh_FEasCP, mesh_NcpElems, mesh_maxNips, mesh_element - use lattice, only: lattice_init - use material - use constitutive, only: constitutive_init,constitutive_state_old,constitutive_state_new + use mesh, only: mesh_init,& + mesh_FEasCP,mesh_element,mesh_NcpElems,mesh_maxNips,FE_Nips + use lattice, only: lattice_init + use material, only: material_init, homogenization_maxNgrains + use constitutive, only: constitutive_init,& + constitutive_state0,constitutive_state + use crystallite + use homogenization implicit none -! + 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,3,3) :: H_bar, H_bar_sym + real(pReal), dimension (3,3) :: ffn,ffn1,Kirchhoff + real(pReal), dimension (3,3,3,3) :: H, H_sym real(pReal), dimension(CPFEM_ngens) :: CPFEM_stress real(pReal), dimension(CPFEM_ngens,CPFEM_ngens) :: CPFEM_jaco real(pReal) Temperature,CPFEM_dt,J_inverse @@ -179,627 +106,116 @@ ! 4: recycling of former results (MARC speciality)& ! 5: record tangent from former converged inc& ! 6: restore tangent from former converged inc - logical CPFEM_updateJaco -! + integer(pInt) e + logical CPFEM_updateJaco + if (.not. CPFEM_init_done) then ! initialization step (three dimensional stress state check missing?) call math_init() call FE_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 material_init() call constitutive_init() - write (6,*) 'call CPFEM init' - call CPFEM_init(Temperature) + call crystallite_init() + call homogenization_init() + call CPFEM_init() CPFEM_init_done = .true. endif -! - if ((.not. parallelExecution) .and. (CPFEM_mode == 3)) CPFEM_mode = 2 -! + cp_en = mesh_FEasCP('elem',CPFEM_en) 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)') & 'theTime',theTime,'theInc',theInc,'theCycle',theCycle,'theLovl',theLovl,& 'mode',CPFEM_mode + write(6,*) '#####################################' endif -! + select case (CPFEM_mode) case (1,2) ! regular computation (with aging of results if mode == 1) if (CPFEM_mode == 1) then ! age results at start of new increment - CPFEM_Lp_old = CPFEM_Lp_new - CPFEM_Fp_old = CPFEM_Fp_new + crystallite_F0 = crystallite_partionedF ! crystallite deformation (_subF is perturbed...) + crystallite_Fp0 = crystallite_Fp ! crystallite plastic deformation + crystallite_Lp0 = crystallite_Lp ! crystallite plastic velocity forall (i = 1:homogenization_maxNgrains,& j = 1:mesh_maxNips, & k = 1:mesh_NcpElems) & - constitutive_state_old(i,j,k)%p = constitutive_state_new(i,j,k)%p - write (6,*) 'results aged.' + constitutive_state0(i,j,k)%p = constitutive_state(i,j,k)%p ! microstructure of crystallites + 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 - if (outdatedFFN1 .or. any(abs(ffn1 - CPFEM_ffn1_bar(:,:,CPFEM_in,cp_en)) > relevantStrain)& - .and. parallelExecution) then - 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) + if (outdatedFFN1 .or. any(abs(ffn1 - materialpoint_F(:,:,CPFEM_in,cp_en)) > relevantStrain)) 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) outdatedFFN1 = .true. - CPFEM_stress_bar(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_cs(1:CPFEM_ngens,CPFEM_in,cp_en) = CPFEM_odd_stress + CPFEM_dcsde(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = CPFEM_odd_jacobian*math_identity2nd(CPFEM_ngens) else if (.not. parallelExecution) then - CPFEM_execution_elem(1) = cp_en - CPFEM_execution_elem(2) = cp_en - CPFEM_execution_IP(1,cp_en) = CPFEM_in - CPFEM_execution_IP(2,cp_en) = CPFEM_in - CPFEM_Temperature(CPFEM_in,cp_en) = Temperature - CPFEM_ffn_bar(:,:,CPFEM_in,cp_en) = ffn - CPFEM_ffn1_bar(:,:,CPFEM_in,cp_en) = ffn1 - call CPFEM_MaterialPoint(CPFEM_updateJaco, CPFEM_dt) + FEsolving_execElem(1) = cp_en + FEsolving_execElem(2) = cp_en + FEsolving_execIP(1,cp_en) = CPFEM_in + FEsolving_execIP(2,cp_en) = CPFEM_in + call materialpoint_stressAndItsTangent(CPFEM_updateJaco, CPFEM_dt) + call materialpoint_postResults(CPFEM_dt) 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. endif ! 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))) - J_inverse = 1.0_pReal/math_det3x3(CPFEM_ffn1_bar(:,:,CPFEM_in, cp_en)) - CPFEM_stress_bar(1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel33to6(J_inverse*Kirchhoff_bar) -! - H_bar = 0.0_pReal + Kirchhoff = math_mul33x33(materialpoint_P(:,:,CPFEM_in, cp_en),transpose(materialpoint_F(:,:,CPFEM_in, cp_en))) + J_inverse = 1.0_pReal/math_det3x3(materialpoint_F(:,:,CPFEM_in, cp_en)) + CPFEM_cs(1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel33to6(J_inverse*Kirchhoff) + + H = 0.0_pReal 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) + & - CPFEM_ffn1_bar(j,m,CPFEM_in,cp_en) * & - CPFEM_ffn1_bar(l,n,CPFEM_in,cp_en) * & - CPFEM_dPdF_bar(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) + & - 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)) + H(i,j,k,l) = H(i,j,k,l) + & + materialpoint_F(j,m,CPFEM_in,cp_en) * & + materialpoint_F(l,n,CPFEM_in,cp_en) * & + materialpoint_dPdF(i,m,k,n,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(j,l) + math_I3(j,l)*Kirchhoff(i,k) + & + 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) & - 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) + H_sym(i,j,k,l)= 0.25_pReal*(H(i,j,k,l)+H(j,i,k,l)+H(i,j,l,k)+H(j,i,l,k)) ! where to use the symmetric version?? + CPFEM_dcsde(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel3333to66(J_inverse*H) endif case (3) ! collect and return odd result - CPFEM_Temperature(CPFEM_in,cp_en) = Temperature - CPFEM_ffn_bar(:,:,CPFEM_in,cp_en) = ffn - CPFEM_ffn1_bar(:,:,CPFEM_in,cp_en) = ffn1 - CPFEM_stress_bar(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) + materialpoint_Temperature(CPFEM_in,cp_en) = Temperature + materialpoint_F0(:,:,CPFEM_in,cp_en) = ffn + materialpoint_F(:,:,CPFEM_in,cp_en) = ffn1 + CPFEM_cs(1:CPFEM_ngens,CPFEM_in,cp_en) = CPFEM_odd_stress + CPFEM_dcsde(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = CPFEM_odd_jacobian*math_identity2nd(CPFEM_ngens) CPFEM_calc_done = .false. 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) - CPFEM_jaco_knownGood = CPFEM_jaco_bar - case (6) ! restore consistent tangent after cutback - CPFEM_jaco_bar = CPFEM_jaco_knownGood + case (5) ! record consistent tangent at beginning of new FE increment (while recycling) + CPFEM_dcsde_knownGood = CPFEM_dcsde + case (6) ! restore consistent tangent after FE cutback + CPFEM_dcsde = CPFEM_dcsde_knownGood end select -! + ! 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_jaco(1:CPFEM_ngens,1:CPFEM_ngens) = CPFEM_jaco_bar(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 + CPFEM_stress(1:CPFEM_ngens) = CPFEM_cs(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 +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 - - 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 !############################################################## \ No newline at end of file diff --git a/trunk/FEsolving.f90 b/trunk/FEsolving.f90 index 577c795b0..aa5c32481 100644 --- a/trunk/FEsolving.f90 +++ b/trunk/FEsolving.f90 @@ -12,6 +12,8 @@ logical :: lastIncConverged = .false.,outdatedByNewInc = .false.,outdatedFFN1 = .false. logical :: symmetricSolver = .false. logical :: parallelExecution = .true. + integer(pInt), dimension(:,:), allocatable :: FEsolving_execIP + integer(pInt), dimension(2) :: FEsolving_execElem CONTAINS @@ -26,7 +28,8 @@ implicit none 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 if (IO_open_inputFile(fileunit)) then @@ -34,11 +37,11 @@ rewind(fileunit) do read (fileunit,'(a1024)',END=100) line - pos = IO_stringPos(line,1) - if( IO_lc(IO_stringValue(line,pos,1)) == 'solver' ) then + positions = IO_stringPos(line,1) + if( IO_lc(IO_stringValue(line,positions,1)) == 'solver' ) then read (fileunit,'(a1024)',END=100) line ! Garbage line - pos = IO_stringPos(line,2) - symmetricSolver = (IO_intValue(line,pos,2) /= 1_pInt) + positions = IO_stringPos(line,2) + symmetricSolver = (IO_intValue(line,positions,2) /= 1_pInt) exit endif enddo diff --git a/trunk/IO.f90 b/trunk/IO.f90 index 57c5d548f..d3c5e5afc 100644 --- a/trunk/IO.f90 +++ b/trunk/IO.f90 @@ -879,6 +879,8 @@ END FUNCTION select case (ID) case (650) msg = 'Polar decomposition failed' + case (600) + msg = 'Crystallite responds elastically' case default msg = 'Unknown warning number...' end select diff --git a/trunk/constitutive.f90 b/trunk/constitutive.f90 index 1616822ec..47e5776c2 100644 --- a/trunk/constitutive.f90 +++ b/trunk/constitutive.f90 @@ -13,11 +13,13 @@ MODULE constitutive use prec implicit none - type(p_vec), dimension(:,:,:), allocatable :: constitutive_state_old, & ! pointer array to old state variables of each grain - constitutive_state_new ! pointer array to new state variables of each grain - integer(pInt), dimension(:,:,:), allocatable :: constitutive_sizeDotState, & ! size of dotState array - constitutive_sizeState, & ! size of state array per grain - constitutive_sizePostResults ! size of postResults array per grain + type(p_vec), dimension(:,:,:), allocatable :: constitutive_state0, & ! pointer array to microstructure at start of FE inc + constitutive_partionedState0, & ! pointer array to microstructure at start of homogenization inc + constitutive_subState0, & ! pointer array to microstructure at start of crystallite inc + constitutive_state ! pointer array to current microstructure (end of converged time step) + integer(pInt), dimension(:,:,:), allocatable :: constitutive_sizeDotState, & ! size of dotState array + constitutive_sizeState, & ! size of state array per grain + constitutive_sizePostResults ! size of postResults array per grain integer(pInt) constitutive_maxSizeDotState,constitutive_maxSizeState,constitutive_maxSizePostResults CONTAINS @@ -44,7 +46,7 @@ subroutine constitutive_init() use constitutive_dislobased 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 @@ -54,44 +56,51 @@ subroutine constitutive_init() close(fileunit) - allocate(constitutive_state_old(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) - allocate(constitutive_state_new(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) + allocate(constitutive_state0(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_sizeState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_sizeState = 0_pInt allocate(constitutive_sizePostResults(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_sizePostResults = 0_pInt do e = 1,mesh_NcpElems ! loop over elements + myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = 1,FE_Nips(mesh_element(2,e)) ! loop over IPs - 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)) select case(phase_constitution(material_phase(g,i,e))) case (constitutive_phenomenological_label) - allocate(constitutive_state_old(g,i,e)%p(constitutive_phenomenological_sizeState(myInstance))) - allocate(constitutive_state_new(g,i,e)%p(constitutive_phenomenological_sizeState(myInstance))) - constitutive_state_new(g,i,e)%p = constitutive_phenomenological_stateInit(myInstance) - constitutive_state_old(g,i,e)%p = constitutive_phenomenological_stateInit(myInstance) + allocate(constitutive_state0(g,i,e)%p(constitutive_phenomenological_sizeState(myInstance))) + allocate(constitutive_partionedState0(g,i,e)%p(constitutive_phenomenological_sizeState(myInstance))) + allocate(constitutive_subState0(g,i,e)%p(constitutive_phenomenological_sizeState(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_sizeState(g,i,e) = constitutive_phenomenological_sizeState(myInstance) constitutive_sizePostResults(g,i,e) = constitutive_phenomenological_sizePostResults(myInstance) case (constitutive_j2_label) - allocate(constitutive_state_old(g,i,e)%p(constitutive_j2_sizeState(myInstance))) - allocate(constitutive_state_new(g,i,e)%p(constitutive_j2_sizeState(myInstance))) - constitutive_state_new(g,i,e)%p = constitutive_j2_stateInit(myInstance) - constitutive_state_old(g,i,e)%p = constitutive_j2_stateInit(myInstance) + allocate(constitutive_state0(g,i,e)%p(constitutive_j2_sizeState(myInstance))) + allocate(constitutive_partionedState0(g,i,e)%p(constitutive_j2_sizeState(myInstance))) + allocate(constitutive_subState0(g,i,e)%p(constitutive_j2_sizeState(myInstance))) + allocate(constitutive_state(g,i,e)%p(constitutive_j2_sizeState(myInstance))) + constitutive_state0(g,i,e)%p = constitutive_j2_stateInit(myInstance) constitutive_sizeDotState(g,i,e) = constitutive_j2_sizeDotState(myInstance) constitutive_sizeState(g,i,e) = constitutive_j2_sizeState(myInstance) constitutive_sizePostResults(g,i,e) = constitutive_j2_sizePostResults(myInstance) case (constitutive_dislobased_label) - allocate(constitutive_state_old(g,i,e)%p(constitutive_dislobased_sizeState(myInstance))) - allocate(constitutive_state_new(g,i,e)%p(constitutive_dislobased_sizeState(myInstance))) - constitutive_state_new(g,i,e)%p = constitutive_dislobased_stateInit(myInstance) - constitutive_state_old(g,i,e)%p = constitutive_dislobased_stateInit(myInstance) + allocate(constitutive_state0(g,i,e)%p(constitutive_dislobased_sizeState(myInstance))) + allocate(constitutive_partionedState0(g,i,e)%p(constitutive_dislobased_sizeState(myInstance))) + allocate(constitutive_subState0(g,i,e)%p(constitutive_dislobased_sizeState(myInstance))) + allocate(constitutive_state(g,i,e)%p(constitutive_dislobased_sizeState(myInstance))) + constitutive_state0(g,i,e)%p = constitutive_dislobased_stateInit(myInstance) constitutive_sizeDotState(g,i,e) = constitutive_dislobased_sizeDotState(myInstance) constitutive_sizeState(g,i,e) = constitutive_dislobased_sizeState(myInstance) constitutive_sizePostResults(g,i,e) = constitutive_dislobased_sizePostResults(myInstance) case default call IO_error(200,material_phase(g,i,e)) ! unknown constitution end select + constitutive_partionedState0(g,i,e)%p = constitutive_state0(g,i,e)%p enddo enddo enddo @@ -99,6 +108,20 @@ subroutine constitutive_init() constitutive_maxSizeState = maxval(constitutive_sizeState) 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 end subroutine @@ -126,11 +149,11 @@ function constitutive_homogenizedC(ipc,ip,el) select case (phase_constitution(material_phase(ipc,ip,el))) 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) - 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) - constitutive_homogenizedC = constitutive_dislobased_homogenizedC(constitutive_state_new,ipc,ip,el) + constitutive_homogenizedC = constitutive_dislobased_homogenizedC(constitutive_state,ipc,ip,el) end select @@ -161,11 +184,11 @@ real(pReal) Temperature select case (phase_constitution(material_phase(ipc,ip,el))) 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) - 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) - call constitutive_dislobased_microstructure(Temperature,constitutive_state_new,ipc,ip,el) + call constitutive_dislobased_microstructure(Temperature,constitutive_state,ipc,ip,el) 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))) 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) - 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) - 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 @@ -241,11 +264,11 @@ function constitutive_dotState(Tstar_v,Temperature,ipc,ip,el) select case (phase_constitution(material_phase(ipc,ip,el))) 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) - 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) - 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 return @@ -278,11 +301,11 @@ pure function constitutive_postResults(Tstar_v,Temperature,dt,ipc,ip,el) constitutive_postResults = 0.0_pReal select case (phase_constitution(material_phase(ipc,ip,el))) case (constitutive_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) - 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) - 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 diff --git a/trunk/crystallite.f90 b/trunk/crystallite.f90 new file mode 100644 index 000000000..513227b38 --- /dev/null +++ b/trunk/crystallite.f90 @@ -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 +!############################################################## \ No newline at end of file diff --git a/trunk/debug.f90 b/trunk/debug.f90 index 4e3b1afb2..de8092789 100644 --- a/trunk/debug.f90 +++ b/trunk/debug.f90 @@ -6,9 +6,10 @@ implicit none - integer(pInt), dimension(nCutback+1) :: debug_cutbackDistribution = 0_pInt - integer(pInt), dimension(nInner) :: debug_InnerLoopDistribution = 0_pInt - integer(pInt), dimension(nOuter) :: debug_OuterLoopDistribution = 0_pInt + integer(pInt), dimension(nLp) :: debug_LpLoopDistribution = 0_pInt + integer(pInt), dimension(nStress) :: debug_StressLoopDistribution = 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_cumDotStateTicks = 0_pInt integer(pInt) :: debug_cumLpCalls = 0_pInt @@ -18,6 +19,24 @@ 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 @@ -47,34 +66,50 @@ dble(debug_cumDotStateTicks)/tickrate/1.0e-6_pReal/debug_cumDotStateCalls write(6,'(a33,x,i12)') 'total CPU ticks :',debug_cumDotStateTicks 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 - write(6,*) 'distribution_InnerLoop :' - do i=1,nInner - if (debug_InnerLoopDistribution(i) /= 0) then - integral = integral + i*debug_InnerLoopDistribution(i) - write(6,*) i,debug_InnerLoopDistribution(i) + write(6,*) + write(6,*) 'distribution_LpLoop :' + do i=1,nLp + if (debug_LpLoopDistribution(i) /= 0) then + integral = integral + i*debug_LpLoopDistribution(i) + write(6,*) i,debug_LpLoopDistribution(i) endif enddo - write(6,*) 'total',sum(debug_InnerLoopDistribution),integral - write(6,*) + write(6,*) 'total',sum(debug_LpLoopDistribution),integral integral = 0_pInt - write(6,*) 'distribution_OuterLoop :' - do i=1,nOuter - if (debug_OuterLoopDistribution(i) /= 0) then - integral = integral + i*debug_OuterLoopDistribution(i) - write(6,*) i,debug_OuterLoopDistribution(i) + write(6,*) + write(6,*) 'distribution_StressLoop :' + do i=1,nStress + if (debug_StressLoopDistribution(i) /= 0) then + integral = integral + i*debug_StressLoopDistribution(i) + write(6,*) i,debug_StressLoopDistribution(i) endif 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,*) END SUBROUTINE diff --git a/trunk/homogenization.f90 b/trunk/homogenization.f90 new file mode 100644 index 000000000..c31daa892 --- /dev/null +++ b/trunk/homogenization.f90 @@ -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 \ No newline at end of file diff --git a/trunk/homogenization_isostrain.f90 b/trunk/homogenization_isostrain.f90 new file mode 100644 index 000000000..38bd0be9a --- /dev/null +++ b/trunk/homogenization_isostrain.f90 @@ -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 + 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 diff --git a/trunk/lattice.f90 b/trunk/lattice.f90 index 554069bfa..9173f3e14 100644 --- a/trunk/lattice.f90 +++ b/trunk/lattice.f90 @@ -566,7 +566,6 @@ subroutine lattice_init() 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 - write(6,*) 'lattice Nstructure',lattice_Nstructure end subroutine @@ -595,7 +594,6 @@ function lattice_initializeStructure(struct,CoverA) integer(pInt) lattice_initializeStructure - write(6,*) 'initialize structure', struct select case(struct(1:3)) ! check first three chars of structure name case ('fcc') myStructure = 1_pInt @@ -705,7 +703,6 @@ function lattice_initializeStructure(struct,CoverA) endif lattice_initializeStructure = myStructure - write(6,*) 'lattice_initializeStructure', myStructure end function diff --git a/trunk/material.f90 b/trunk/material.f90 index 8f76b4837..2be57113d 100644 --- a/trunk/material.f90 +++ b/trunk/material.f90 @@ -51,8 +51,8 @@ integer(pInt), dimension(:), allocatable :: homogenization_Ngrains, & integer(pInt), dimension(:,:), allocatable :: microstructure_phase, & ! phase 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 :: material_volFrac ! vol fraction of grain within phase (?) -integer(pInt), dimension(:,:,:), allocatable :: material_phase ! phase of each grain,IP,element +integer(pInt), dimension(:,:,:), allocatable :: material_volFrac, & ! vol fraction of grain within phase (?) + 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 :: texture_Gauss, & ! data of each Gauss component texture_Fiber ! data of each Fiber component @@ -74,13 +74,9 @@ subroutine material_init() integer(pInt) i 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) - write(6,*) 'parsing microstrcuture...' call material_parseMicrostructure(fileunit,material_partMicrostructure) - write(6,*) 'parsing texture...' call material_parseTexture(fileunit,material_partTexture) - write(6,*) 'parsing phase...' call material_parsePhase(fileunit,material_partPhase) close(fileunit) @@ -104,9 +100,7 @@ subroutine material_init() write (6,'(a32,x,i4)') microstructure_name(i),microstructure_Nconstituents(i) enddo - write(6,*) 'populating grains...' call material_populateGrains() - write(6,*) 'populating grains finished...' end subroutine @@ -127,14 +121,19 @@ subroutine material_parseHomogenization(file,myPart) character(len=64) tag character(len=1024) line - Nsections= IO_countSections(file,myPart) + Nsections = IO_countSections(file,myPart) material_Nhomogenization = Nsections + + write (6,*) 'homogenization sections found',material_Nhomogenization allocate(homogenization_name(Nsections)); homogenization_name = '' allocate(homogenization_type(Nsections)); homogenization_type = '' allocate(homogenization_typeInstance(Nsections)); homogenization_typeInstance = 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) + write(6,*) 'count of (output)',homogenization_Noutput rewind(file) line = '' @@ -434,6 +433,7 @@ subroutine material_populateGrains() integer(pInt), dimension (:,:), allocatable :: Ngrains integer(pInt), dimension (microstructure_maxNconstituents) :: NgrainsOfConstituent + real(pReal), dimension (:,:), allocatable :: volume real(pReal), dimension (:), allocatable :: volFracOfGrain, phaseOfGrain real(pReal), dimension (:,:), allocatable :: orientationOfGrain 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(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 homog = mesh_element(3,e) micro = mesh_element(4,e) @@ -458,6 +459,7 @@ subroutine material_populateGrains() if (micro < 1 .or. micro > material_Nmicrostructure) & ! out of bounds call IO_error(140,e,0,0) 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 allocate(volFracOfGrain(maxval(Ngrains))) ! reserve memory for maximum case @@ -480,16 +482,15 @@ subroutine material_populateGrains() 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 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 endif enddo - write (6,*) 'now at grain count',grain ! ---------------------------------------------------------------------------- NgrainsOfConstituent = 0_pInt forall (i = 1:microstructure_Nconstituents(micro)) & NgrainsOfConstituent(i) = nint(microstructure_fraction(i,micro) * myNgrains, pInt) - write (6,*) 'NgrainsOfConstituent',NgrainsOfConstituent do while (sum(NgrainsOfConstituent) /= myNgrains) ! total grain count over constituents wrong? sgn = sign(1_pInt, myNgrains - sum(NgrainsOfConstituent)) ! direction of required change extreme = 0.0_pReal @@ -502,8 +503,6 @@ subroutine material_populateGrains() enddo NgrainsOfConstituent(t) = NgrainsOfConstituent(t) + sgn ! change that by one end do - write (6,*) 'fixed NgrainsOfConstituent',NgrainsOfConstituent - ! ---------------------------------------------------------------------------- phaseOfGrain = 0_pInt orientationOfGrain = 0.0_pReal @@ -515,26 +514,21 @@ subroutine material_populateGrains() 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) - 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 ! --------- if (texture_ODFfile(textureID) == '') then ! dealing with texture 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 orientationOfGrain(:,grain+constituentGrain+g) = & math_sampleGaussOri(texture_Gauss(1:3,t,textureID),& texture_Gauss( 4,t,textureID)) enddo constituentGrain = constituentGrain + int(myNorientations*texture_Gauss(5,t,textureID)) - write (6,*) 'now at constituent grain',constituentGrain enddo 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 orientationOfGrain(:,grain+constituentGrain+g) = & math_sampleFiberOri(texture_Fiber(1:2,t,textureID),& @@ -542,14 +536,11 @@ subroutine material_populateGrains() texture_Fiber( 5,t,textureID)) enddo constituentGrain = constituentGrain + int(myNorientations*texture_fiber(6,t,textureID)) - write (6,*) 'now at constituent grain',constituentGrain enddo - write (6,*) 'looping',constituentGrain+1,myNorientations do j = constituentGrain+1,myNorientations ! fill remainder with random orientationOfGrain(:,grain+j) = math_sampleRandomOri() enddo - write (6,*) 'done...' ! --------- else ! hybrid IA ! --------- @@ -592,9 +583,6 @@ subroutine material_populateGrains() !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 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 @@ -603,13 +591,9 @@ subroutine material_populateGrains() material_phase(g,i,e) = phaseOfGrain(grain+(i-1)*dGrains+g) material_EulerAngles(:,g,i,e) = orientationOfGrain(:,grain+(i-1)*dGrains+g) end forall - !do i = 1,FE_Nips(mesh_element(2,e)) - ! write(6,'(3(i10,x),e10.3,x,i10,x,3(f10.1,x))') e, i, dGrains, sum(material_volFrac(:,i,e)), & - ! sum(material_phase(:,i,e)), & - ! sum(material_EulerAngles(1,:,i,e)), & - ! sum(material_EulerAngles(2,:,i,e)), & - ! sum(material_EulerAngles(3,:,i,e)) - !end do + write (6,*) e + write (6,*) material_phase(:,:,e) + write (6,*) material_EulerAngles(:,:,:,e) grain = grain + FE_Nips(mesh_element(2,e)) * dGrains ! wind forward by Nips*NgrainsPerIP endif enddo diff --git a/trunk/math.f90 b/trunk/math.f90 index 2791f2f8d..39f20a994 100644 --- a/trunk/math.f90 +++ b/trunk/math.f90 @@ -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 error = .false. -!!$OMP CRITICAL (evilmatmul) - - ce=math_mul33x33(transpose(fe),fe) -!!$OMP END CRITICAL (evilmatmul) + ce = math_mul33x33(transpose(FE),FE) CALL math_spectral1(CE,EW1,EW2,EW3,EB1,EB2,EB3) U=DSQRT(EW1)*EB1+DSQRT(EW2)*EB2+DSQRT(EW3)*EB3 call math_invert3x3(U,UI,det,error) - if (.not. error) R = math_mul33x33(fe,ui) + if (.not. error) R = math_mul33x33(FE,UI) return diff --git a/trunk/mesh.f90 b/trunk/mesh.f90 index adf8a19ea..64e5a3366 100644 --- a/trunk/mesh.f90 +++ b/trunk/mesh.f90 @@ -1009,7 +1009,6 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements END SUBROUTINE - !******************************************************************** ! get count of elements, nodes, and cp elements in mesh ! for subsequent array allocations @@ -1685,18 +1684,18 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit) !$OMP CRITICAL (write2out) - !write (6,*) - !write (6,*) "Input Parser: IP NEIGHBORHOOD" - !write (6,*) - !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 - ! t = mesh_element(2,e) ! get elemType - ! do i = 1,FE_Nips(t) ! loop over IPs of elem - ! 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) - ! enddo - ! enddo - !enddo + write (6,*) + write (6,*) "Input Parser: IP NEIGHBORHOOD" + write (6,*) + 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 + t = mesh_element(2,e) ! get elemType + do i = 1,FE_Nips(t) ! loop over IPs of elem + 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) + enddo + enddo + enddo write (6,*) write (6,*) "Input Parser: ELEMENT VOLUME" write (6,*) @@ -1706,9 +1705,9 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit) do e = 1,mesh_NcpElems do i = 1,FE_Nips(mesh_element(2,e)) write (6,"(i5,x,i5,x,e15.8)") e,i,mesh_IPvolume(i,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) - ! end do + 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) + end do end do end do !write (6,*) @@ -1746,11 +1745,11 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit) write (6,*) mesh_maxValStateVar(1), " : maximum homogenization index" write (6,*) mesh_maxValStateVar(2), " : maximum microstructure index" write (6,*) - write (fmt,"(a,i5,a)") "(9(x),a1,x,",mesh_maxValStateVar(2),"(i8))" - 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)") "(9(x),a2,x,",mesh_maxValStateVar(2),"(i8))" + write (6,fmt) "+-",math_range(mesh_maxValStateVar(2)) + 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 - 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 write (6,*) !$OMP END CRITICAL (write2out) diff --git a/trunk/mpie_cpfem_marc.f90 b/trunk/mpie_cpfem_marc.f90 index 642982929..e846b71ba 100644 --- a/trunk/mpie_cpfem_marc.f90 +++ b/trunk/mpie_cpfem_marc.f90 @@ -37,19 +37,18 @@ 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 "mesh.f90" ! uses prec, math, IO, FEsolving 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_j2.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 "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 - 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 @@ -75,51 +74,6 @@ ! 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 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) @@ -131,12 +85,49 @@ !2 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 FEsolving use CPFEM, only: CPFEM_general use math, only: invnrmMandel use debug -! implicit none ! ** Start of generated type statements ** @@ -146,7 +137,7 @@ integer(pInt) ndi, ndm, ngens, nn, nnode, nshear real(pReal) s, strechn, strechn1, t ! ** 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,*),& 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 -! write(6,'(3(3(f10.3,x),/))') ffn1(:,1),ffn1(:,2),ffn1(:,3) - if (inc == 0) then cycleCounter = 4 else @@ -169,13 +158,7 @@ outdatedFFN1 = .false. write (6,*) n(1),nn,'cycleCounter',cycleCounter call debug_info() ! output of debugging/performance statistics of former - debug_cutbackDistribution = 0_pInt ! initialize debugging data - debug_InnerLoopDistribution = 0_pInt - debug_OuterLoopDistribution = 0_pInt - debug_cumLpTicks = 0 - debug_cumDotStateTicks = 0 - debug_cumLpCalls = 0_pInt - debug_cumDotStateCalls = 0_pInt + call debug_reset() endif endif if (cptim > theTime .or. theInc /= inc) then ! reached convergence @@ -213,48 +196,47 @@ 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 !******************************************************************** ! ! 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 CPFEM, only: CPFEM_results, CPFEM_Nresults - use constitutive, only: constitutive_maxSizePostResults use mesh, only: mesh_FEasCP + use homogenization, only: materialpoint_results 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_maxSizePostResults)+1_pInt,& - (jpltcd-1_pInt)/(CPFEM_Nresults+constitutive_maxSizePostResults)+1_pInt,& - nn, mesh_FEasCP('elem', m)) + + v = materialpoint_results(jpltcd,nn,mesh_FEasCP('elem', m)) return - END SUBROUTINE -! -! + +end subroutine + + + ! subroutine utimestep(timestep,timestepold,icall,time,timeloadcase) !******************************************************************** ! This routine modifies the addaptive time step of Marc diff --git a/trunk/prec.f90 b/trunk/prec.f90 index 5aef6e736..7a617f01d 100644 --- a/trunk/prec.f90 +++ b/trunk/prec.f90 @@ -20,14 +20,15 @@ ! *** Numerical parameters *** 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 - integer(pInt), parameter :: nOuter = 20_pInt ! outer loop limit 20 - integer(pInt), parameter :: nInner = 200_pInt ! inner loop limit 200 - real(pReal), parameter :: reltol_Outer = 1.0e-5_pReal ! relative tolerance in outer loop (state) - real(pReal), parameter :: reltol_Inner = 1.0e-6_pReal ! relative tolerance in inner loop (Lp) - real(pReal), parameter :: abstol_Inner = 1.0e-8_pReal ! absolute tolerance in inner loop (Lp) + integer(pInt), parameter :: nReg = 1_pInt ! regularization attempts for Jacobi inversion + integer(pInt), parameter :: nHomog = 10_pInt ! homogenization loop limit + integer(pInt), parameter :: nCryst = 10_pInt ! crystallite loop limit (state update) + integer(pInt), parameter :: nStress = 20_pInt ! stress loop limit + 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 :: resAbsol = 1.0e+2_pReal ! absolute tolerance of residual in GIA iteration (corresponds to ~1 Pa)