diff --git a/trunk/CPFEM_Taylor.f90 b/trunk/CPFEM_Taylor.f90 index 1af6a5326..2f9b010df 100644 --- a/trunk/CPFEM_Taylor.f90 +++ b/trunk/CPFEM_Taylor.f90 @@ -1,157 +1,192 @@ -!############################################################## - MODULE CPFEM -!############################################################## -! *** CPFEM engine *** -! - use prec, only: pReal,pInt - implicit none -! -! **************************************************************** -! *** General variables for the material behaviour calculation *** -! **************************************************************** - real(pReal), dimension (:,:), allocatable :: CPFEM_Temperature - real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_ffn_bar - real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_ffn1_bar - real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_PK1_bar - real(pReal), dimension (:,:,:,:,:,:),allocatable :: CPFEM_dPdF_bar - 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 - real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_Fp_old - real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_Fp_new - real(pReal), parameter :: CPFEM_odd_stress = 1e15_pReal, CPFEM_odd_jacobian = 1e50_pReal - integer(pInt) :: CPFEM_Nresults = 4_pInt ! three Euler angles plus volume fraction - logical :: CPFEM_init_done = .false. ! remember if init has been done already - logical :: CPFEM_calc_done = .false. ! remember if first IP has already calced the results -! - CONTAINS -! -!********************************************************* -!*** 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 mesh - use constitutive -! - 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_ffn1_bar (3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_ffn1_bar = CPFEM_ffn_bar - allocate(CPFEM_PK1_bar (3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_PK1_bar = 0.0_pReal - allocate(CPFEM_dPdF_bar(3,3,3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dPdF_bar = 0.0_pReal - allocate(CPFEM_stress_bar(6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_stress_bar = 0.0_pReal - allocate(CPFEM_jaco_bar(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_jaco_bar = 0.0_pReal - allocate(CPFEM_jaco_knownGood(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_jaco_knownGood = 0.0_pReal -! -! *** User defined results !!! MISSING incorporate consti_Nresults *** - allocate(CPFEM_results(CPFEM_Nresults+constitutive_maxNresults,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) - CPFEM_results = 0.0_pReal -! -! *** Plastic velocity gradient *** - allocate(CPFEM_Lp(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Lp = 0.0_pReal - -! *** Plastic deformation gradient at (t=t0) and (t=t1) *** - allocate(CPFEM_Fp_new(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Fp_new = 0.0_pReal - allocate(CPFEM_Fp_old(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) - forall (e=1:mesh_NcpElems,i=1:mesh_maxNips,g=1:constitutive_maxNgrains) & - CPFEM_Fp_old(:,:,g,i,e) = math_EulerToR(constitutive_EulerAngles(:,g,i,e)) ! plastic def gradient reflects init orientation -! -! *** Output to MARC output file *** -!$OMP CRITICAL (write2out) - write(6,*) - write(6,*) 'CPFEM Initialization' - write(6,*) - write(6,*) 'CPFEM_Temperature: ', shape(CPFEM_Temperature) - write(6,*) 'CPFEM_ffn_bar: ', shape(CPFEM_ffn_bar) - write(6,*) 'CPFEM_ffn1_bar: ', shape(CPFEM_ffn1_bar) - write(6,*) 'CPFEM_PK1_bar: ', shape(CPFEM_PK1_bar) - write(6,*) 'CPFEM_dPdF_bar: ', shape(CPFEM_dPdF_bar) - 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: ', shape(CPFEM_Lp) - write(6,*) 'CPFEM_Fp_old: ', shape(CPFEM_Fp_old) - write(6,*) 'CPFEM_Fp_new: ', shape(CPFEM_Fp_new) - write(6,*) - call flush(6) -!$OMP END CRITICAL (write2out) - return -! - END SUBROUTINE -! -! -!*********************************************************************** -!*** perform initialization at first call, update variables and *** -!*** call the actual material model *** -! -! CPFEM_mode computation mode (regular, collection, recycle) -! ffn deformation gradient for t=t0 -! ffn1 deformation gradient for t=t1 -! Temperature temperature -! CPFEM_dt time increment -! CPFEM_en element number -! CPFEM_in intergration point number -! CPFEM_stress stress vector in Mandel notation -! CPFEM_updateJaco flag to initiate computation of Jacobian -! 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) -! note: CPFEM_stress = Cauchy stress cs(6) and CPFEM_jaco = Consistent tangent dcs/de -! - use prec, only: pReal,pInt - use FEsolving - use debug - use math - use mesh, only: mesh_init,mesh_FEasCP, mesh_NcpElems, FE_Nips, FE_mapElemtype, mesh_element - use lattice, only: lattice_init - use constitutive, only: constitutive_init,constitutive_state_old,constitutive_state_new,material_Cslip_66 - implicit none -! - integer(pInt) CPFEM_en, CPFEM_in, cp_en, CPFEM_ngens, i,j,k,l,m,n, e - real(pReal), dimension (3,3) :: ffn,ffn1,Kirchhoff_bar - real(pReal), dimension (3,3,3,3) :: H_bar, H_bar_sym - real(pReal), dimension(CPFEM_ngens) :: CPFEM_stress - real(pReal), dimension(CPFEM_ngens,CPFEM_ngens) :: CPFEM_jaco, odd_jaco - real(pReal) Temperature,CPFEM_dt,J_inverse - integer(pInt) CPFEM_mode ! 1: regular computation with aged results& - ! 2: regular computation& - ! 3: collection of FEM data& - ! 4: recycling of former results (MARC speciality)& - ! 5: record tangent from former converged inc& - ! 6: restore tangent from former converged inc - logical CPFEM_updateJaco -! - if (.not. CPFEM_init_done) then ! initialization step (three dimensional stress state check missing?) - call math_init() - call mesh_init() - call lattice_init() - call constitutive_init() - call CPFEM_init(Temperature) - CPFEM_init_done = .true. - endif -! - cp_en = mesh_FEasCP('elem',CPFEM_en) - if (cp_en == 1 .and. CPFEM_in == 1) & - 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 -! - select case (CPFEM_mode) +!############################################################## + MODULE CPFEM +!############################################################## +! *** CPFEM engine *** +! + use prec, only: pReal,pInt + implicit none +! +! **************************************************************** +! *** 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_Fe1 + real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_Tstar_v + real(pReal), parameter :: CPFEM_odd_stress = 1e15_pReal, CPFEM_odd_jacobian = 1e50_pReal + integer(pInt) :: CPFEM_Nresults = 4_pInt ! three Euler angles plus volume fraction + 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 +! +! *** Solution at single crystallite level *** +! + logical, dimension (:,:,:),allocatable :: crystallite_converged !individual covergence flag per grain +! + CONTAINS +! +!********************************************************* +!*** 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 mesh + use constitutive +! + 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,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) + forall(g=1:constitutive_maxNgrains,e=1:mesh_NcpElems,i=1:mesh_maxNips) CPFEM_ffn(:,:,g,i,e) = math_I3 + allocate(CPFEM_ffn1_bar(3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_ffn1_bar = CPFEM_ffn_bar + allocate(CPFEM_ffn1(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_ffn1 = CPFEM_ffn + allocate(CPFEM_PK1_bar(3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_PK1_bar = 0.0_pReal + allocate(CPFEM_PK1(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_PK1 = 0.0_pReal + allocate(CPFEM_dPdF_bar(3,3,3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dPdF_bar = 0.0_pReal + allocate(CPFEM_dPdF_bar_old(3,3,3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dPdF_bar_old = 0.0_pReal + allocate(CPFEM_dPdF(3,3,3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dPdF = 0.0_pReal + allocate(CPFEM_stress_bar(6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_stress_bar = 0.0_pReal + allocate(CPFEM_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 !!! MISSING incorporate consti_Nresults *** + allocate(CPFEM_results(CPFEM_Nresults+constitutive_maxNresults,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) + CPFEM_results = 0.0_pReal +! +! *** Plastic velocity gradient *** + allocate(CPFEM_Lp_old(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Lp_old = 0.0_pReal + allocate(CPFEM_Lp_new(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Lp_new = 0.0_pReal + +! *** Plastic deformation gradient at (t=t0) and (t=t1) *** + allocate(CPFEM_Fp_new(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Fp_new = 0.0_pReal + allocate(CPFEM_Fp_old(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) + forall (e=1:mesh_NcpElems,i=1:mesh_maxNips,g=1:constitutive_maxNgrains) & + CPFEM_Fp_old(:,:,g,i,e) = math_EulerToR(constitutive_EulerAngles(:,g,i,e)) ! plastic def gradient reflects init orientation +! *** Elastic deformation gradient at (t=t1) *** + allocate(CPFEM_Fe1(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Fe1 = 0.0_pReal +! *** Stress vector at (t=t1) *** + allocate(CPFEM_Tstar_v(6,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Tstar_v = 0.0_pReal +! +! *** Output to MARC output file *** +!$OMP CRITICAL (write2out) + write(6,*) + write(6,*) 'CPFEM Initialization' + 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_Fe1: ', shape(CPFEM_Fe1) + write(6,*) 'CPFEM_Tstar_v: ', shape(CPFEM_Tstar_v) + write(6,*) + call flush(6) +!$OMP END CRITICAL (write2out) + return +! + END SUBROUTINE +! +! +!*********************************************************************** +!*** perform initialization at first call, update variables and *** +!*** call the actual material model *** +! +! CPFEM_mode computation mode (regular, collection, recycle) +! ffn deformation gradient for t=t0 +! ffn1 deformation gradient for t=t1 +! Temperature temperature +! CPFEM_dt time increment +! CPFEM_en element number +! CPFEM_in intergration point number +! CPFEM_stress stress vector in Mandel notation +! CPFEM_updateJaco flag to initiate computation of Jacobian +! 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) +! note: CPFEM_stress = Cauchy stress cs(6) and CPFEM_jaco = Consistent tangent dcs/de +! + use prec, only: pReal,pInt + use FEsolving + use debug + use math + use mesh, only: mesh_init,mesh_FEasCP, mesh_NcpElems, FE_Nips, FE_mapElemtype, mesh_element + use lattice, only: lattice_init + use constitutive, only: constitutive_init,constitutive_state_old,constitutive_state_new,material_Cslip_66 + 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(CPFEM_ngens) :: CPFEM_stress + real(pReal), dimension(CPFEM_ngens,CPFEM_ngens) :: CPFEM_jaco, odd_jaco + real(pReal) Temperature,CPFEM_dt,J_inverse + integer(pInt) CPFEM_mode ! 1: regular computation with aged results& + ! 2: regular computation& + ! 3: collection of FEM data& + ! 4: recycling of former results (MARC speciality)& + ! 5: record tangent from former converged inc& + ! 6: restore tangent from former converged inc + logical CPFEM_updateJaco +! + if (.not. CPFEM_init_done) then ! initialization step (three dimensional stress state check missing?) + call math_init() + call mesh_init() + call lattice_init() + call constitutive_init() + call crystallite_init() + call CPFEM_init(Temperature) + CPFEM_init_done = .true. + endif +! + cp_en = mesh_FEasCP('elem',CPFEM_en) + if (cp_en == 1 .and. CPFEM_in == 1) then +!$OMP CRITICAL (write2out) + 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 +!$OMP END CRITICAL (write2out) + endif +! + select case (CPFEM_mode) case (2,1) ! regular computation (with aging of results) if (any(abs(ffn1 - CPFEM_ffn1_bar(:,:,CPFEM_in,cp_en)) > relevantStrain)) then CPFEM_stress_bar(1:CPFEM_ngens,:,:) = CPFEM_odd_stress @@ -160,6 +195,7 @@ CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,j,i) = odd_jaco outdatedFFN1 = .true. if (.not. CPFEM_calc_done .AND.CPFEM_mode == 1) then + CPFEM_Lp_old = CPFEM_Lp_new CPFEM_Fp_old = CPFEM_Fp_new constitutive_state_old = constitutive_state_new endif @@ -167,45 +203,49 @@ write(6,*) 'WARNING: FFN1 changed for ip', CPFEM_in, 'of element', cp_en !$OMP END CRITICAL (write2out) return - else - if (.not. CPFEM_calc_done) then ! puuh, me needs doing all the work... - if (CPFEM_mode == 1) then ! age results at start of new increment - CPFEM_Fp_old = CPFEM_Fp_new - constitutive_state_old = constitutive_state_new - endif - debug_cutbackDistribution = 0_pInt ! initialize debugging data - debug_InnerLoopDistribution = 0_pInt - debug_OuterLoopDistribution = 0_pInt -! -!$OMP PARALLEL DO - do e=1,mesh_NcpElems ! ## this shall be done in a parallel loop in the future ## - do i=1,FE_Nips(mesh_element(2,e)) ! iterate over all IPs of this element's type - debugger = (e==1 .and. i==1) ! switch on debugging for first IP in first element - call CPFEM_MaterialPoint(CPFEM_updateJaco, CPFEM_dt, i, e) - enddo - enddo -!$OMP END PARALLEL DO - call debug_info() ! output of debugging/performance statistics + else + if (.not. CPFEM_calc_done) then ! puuh, me needs doing all the work... + if (CPFEM_mode == 1) then ! age results at start of new increment + CPFEM_Lp_old = CPFEM_Lp_new + CPFEM_Fp_old = CPFEM_Fp_new + constitutive_state_old = constitutive_state_new + endif + debug_cutbackDistribution = 0_pInt ! initialize debugging data + debug_InnerLoopDistribution = 0_pInt + debug_OuterLoopDistribution = 0_pInt +! +!****************************************************************************************************** +! parallelization moved down to material point and single crystallite +!****************************************************************************************************** +!!$OMP PARALLEL DO +! do e=1,mesh_NcpElems ! ## this shall be done in a parallel loop in the future ## +! do i=1,FE_Nips(mesh_element(2,e)) ! iterate over all IPs of this element's type +! debugger = (e==1 .and. i==1) ! switch on debugging for first IP in first element + call CPFEM_MaterialPoint(CPFEM_updateJaco, CPFEM_dt, 0, 0) +! enddo +! enddo +!!$OMP END PARALLEL DO + call debug_info() ! output of debugging/performance statistics CPFEM_calc_done = .true. ! now calc is done - endif -! translate from P and dP/dF to CS and dCS/dE + endif +! translate from P and dP/dF to CS and dCS/dE !!$OMP CRITICAL (evilmatmul) - Kirchhoff_bar = math_mul33x33(CPFEM_PK1_bar(:,:,CPFEM_in, cp_en),transpose(CPFEM_ffn1_bar(:,:,CPFEM_in, cp_en))) + Kirchhoff_bar = math_mul33x33(CPFEM_PK1_bar(:,:,CPFEM_in, cp_en),transpose(CPFEM_ffn1_bar(:,:,CPFEM_in, cp_en))) !!$OMP END CRITICAL (evilmatmul) - 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 - 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) + & + 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 + 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)) forall(i=1:3,j=1:3,k=1:3,l=1:3) & - H_bar_sym(i,j,k,l)= 0.25_pReal*(H_bar(i,j,k,l) +H_bar(j,i,k,l) +H_bar(i,j,l,k) +H_bar(j,i,l,k)) + H_bar_sym(i,j,k,l)= 0.25_pReal*(H_bar(i,j,k,l) +H_bar(j,i,k,l) +H_bar(i,j,l,k) +H_bar(j,i,l,k)) CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel3333to66(J_inverse*H_bar) ! if (CPFEM_in==8 .and. cp_en==80) then ! do e=1,80 @@ -217,16 +257,16 @@ ! write(6,*) CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,i,e) ! enddo ! enddo -! endif +! endif ! - 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) - CPFEM_calc_done = .false. + 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) + CPFEM_calc_done = .false. ! if (CPFEM_in==8 .and. cp_en==80) then ! do e=1,80 ! do i=1,8 @@ -236,103 +276,721 @@ ! enddo ! enddo ! endif - - case (4) ! do nothing since we can recycle the former results (MARC specialty) - case (5) ! record consistent tangent at beginning of new increment - CPFEM_jaco_knownGood = CPFEM_jaco_bar - case (6) ! restore consistent tangent after cutback - CPFEM_jaco_bar = CPFEM_jaco_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) + + case (4) ! do nothing since we can recycle the former results (MARC specialty) + case (5) ! record consistent tangent at beginning of new increment + CPFEM_jaco_knownGood = CPFEM_jaco_bar + case (6) ! restore consistent tangent after cutback + CPFEM_jaco_bar = CPFEM_jaco_knownGood + end select ! - return -! - END SUBROUTINE -! -! -!********************************************************** -!*** calculate the material point behaviour *** -!********************************************************** - SUBROUTINE CPFEM_MaterialPoint(& - updateJaco,& ! flag to initiate Jacobian updating - CPFEM_dt,& ! Time increment (dt) - CPFEM_in,& ! Integration point number - cp_en) ! Element number -! - use prec - use FEsolving, only: theCycle - use debug - use math, only: math_pDecomposition,math_RtoEuler,inDeg,math_I3,math_invert3x3,math_permut,math_invert,math_delta - use IO, only: IO_error - use mesh, only: mesh_element - use crystallite - use constitutive - implicit none -! - character(len=128) msg - integer(pInt) cp_en,CPFEM_in,grain - logical updateJaco,error - real(pReal) CPFEM_dt,volfrac - real(pReal), dimension(3,3) :: U,R,Fe1 - real(pReal), dimension(3,3) :: PK1 - real(pReal), dimension(3,3,3,3) :: dPdF,dPdF_bar_old -! - CPFEM_PK1_bar(:,:,CPFEM_in,cp_en) = 0.0_pReal ! zero out average first PK stress - if (updateJaco) then - dPdF_bar_old = CPFEM_dPdF_bar(:,:,:,:,CPFEM_in,cp_en) ! remember former average consistent tangent - CPFEM_dPdF_bar(:,:,:,:,CPFEM_in,cp_en) = 0.0_pReal ! zero out avg consistent tangent for later assembly - endif - - do grain = 1,texture_Ngrains(mesh_element(4,cp_en)) - dPdF = dPdF_bar_old ! preguess consistent tangent of grain with avg - call SingleCrystallite(msg,PK1,dPdF,& - CPFEM_results(CPFEM_Nresults+1:CPFEM_Nresults+constitutive_Nresults(grain,CPFEM_in,cp_en),& - grain,CPFEM_in,cp_en),& - CPFEM_Lp(:,:,grain,CPFEM_in,cp_en),& - CPFEM_Fp_new(:,:,grain,CPFEM_in,cp_en),Fe1,constitutive_state_new(:,grain,CPFEM_in,cp_en),& ! output up to here - CPFEM_dt,cp_en,CPFEM_in,grain,updateJaco,& - CPFEM_Temperature(CPFEM_in,cp_en),& - CPFEM_ffn1_bar(:,:,CPFEM_in,cp_en),CPFEM_ffn_bar(:,:,CPFEM_in,cp_en),& - CPFEM_Fp_old(:,:,grain,CPFEM_in,cp_en),constitutive_state_old(:,grain,CPFEM_in,cp_en)) - - if (msg /= 'ok') then ! solution not reached --> exit +! 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) + CPFEM_in,& ! Integration point number + cp_en) ! Element number +! + use prec + use FEsolving, only: theCycle + use debug + use math, only: math_pDecomposition,math_RtoEuler,inDeg,math_I3,math_invert3x3,math_permut,math_invert,math_delta + use IO, only: IO_error + use mesh, only: mesh_element, mesh_NcpElems, FE_Nips +! use crystallite + use constitutive + implicit none +! + integer(pInt) cp_en,CPFEM_in,g,i,e + integer(pInt) el_start, el_end, ip_start, ip_end + logical updateJaco,error + real(pReal) CPFEM_dt,volfrac + real(pReal), dimension(3,3) :: U,R !,Fe1 +! real(pReal), dimension(3,3) :: PK1 +! real(pReal), dimension(3,3,3,3) :: dPdF,dPdF_bar_old +! + CPFEM_PK1_bar = 0.0_pReal ! zero out average first PK stress +!initialize element loop + if (cp_en /= 0_pInt) then + el_start = cp_en + el_end = cp_en + else + el_start = 1_pInt + el_end = mesh_NcpElems + endif +! prescribe FFN and FFN1 depending on homogenization scheme +!$OMP PARALLEL DO + do e=el_start,el_end + if(CPFEM_in /= 0_pInt) then + ip_start = CPFEM_in + ip_end = CPFEM_in + else + ip_start = 1 + ip_end = FE_Nips(mesh_element(2,e)) ! iterate over all IPs of this element's type + endif + do i=ip_start,ip_end + do g=1,texture_Ngrains(mesh_element(4,e)) + CPFEM_ffn(:,:,g,i,e) = CPFEM_ffn_bar(:,:,i,e) !Taylor homogenization + CPFEM_ffn1(:,:,g,i,e) = CPFEM_ffn1_bar(:,:,i,e) !Taylor homogenization + end do + end do + end do +!$OMP END PARALLEL DO +! calculate stress, update state and update jacobian in case needed for all or one ip + if (updateJaco) then + CPFEM_dPdF_bar_old = CPFEM_dPdF_bar ! remember former average consistent tangent + CPFEM_dPdF_bar = 0.0_pReal ! zero out avg consistent tangent for later assembly + endif + call SingleCrystallite(updateJaco,CPFEM_dt,el_start,el_end,CPFEM_in) +!****************************************************************************************************** +! check convergence of homogenization in case needed +!****************************************************************************************************** +! calculate average quantities per ip and post results +!$OMP PARALLEL DO + do e=el_start,el_end + if(CPFEM_in /= 0_pInt) then + ip_start = CPFEM_in + ip_end = CPFEM_in + else + ip_start = 1 + ip_end = FE_Nips(mesh_element(2,e)) ! iterate over all IPs of this element's type + endif + do i=ip_start,ip_end + do g=1,texture_Ngrains(mesh_element(4,e)) + volfrac = constitutive_matVolFrac(g,i,e)*constitutive_texVolFrac(g,i,e) + CPFEM_PK1_bar(:,:,i,e) = CPFEM_PK1_bar(:,:,i,e) + volfrac * CPFEM_PK1(:,:,g,i,e) + if (updateJaco) CPFEM_dPdF_bar(:,:,:,:,i,e) = & + CPFEM_dPdF_bar(:,:,:,:,i,e) + volfrac * CPFEM_dPdF(:,:,:,:,g,i,e) ! add up crystallite stiffnesses + ! (may have "holes" corresponding + ! to former avg tangent) +! update results plotted in MENTAT + call math_pDecomposition(CPFEM_Fe1(:,:,g,i,e),U,R,error) ! polar decomposition + if (error) then !$OMP CRITICAL (write2out) - write(6,*) 'grain loop failed to converge @ EL:',cp_en,' IP:',CPFEM_in + write(6,*) 'polar decomposition of', CPFEM_Fe1(:,:,g,i,e) + write(6,*) 'Grain: ',g + write(6,*) 'Integration point: ',i + write(6,*) 'Element: ',mesh_element(1,e) !$OMP END CRITICAL (write2out) - call IO_error(600) - return - endif - - volfrac = constitutive_matVolFrac(grain,CPFEM_in,cp_en)*constitutive_texVolFrac(grain,CPFEM_in,cp_en) - CPFEM_PK1_bar(:,:,CPFEM_in,cp_en) = CPFEM_PK1_bar(:,:,CPFEM_in,cp_en) + volfrac*PK1 - if (updateJaco) CPFEM_dPdF_bar(:,:,:,:,CPFEM_in,cp_en) = & - CPFEM_dPdF_bar(:,:,:,:,CPFEM_in,cp_en) + volfrac*dPdF ! add up crystallite stiffnesses - ! (may have "holes" corresponding - ! to former avg tangent) -! -! update results plotted in MENTAT - call math_pDecomposition(Fe1,U,R,error) ! polar decomposition - if (error) then -!$OMP CRITICAL (write2out) - write(6,*) 'polar decomposition of', Fe1 - write(6,*) 'Grain: ',grain - write(6,*) 'Integration point: ',CPFEM_in - write(6,*) 'Element: ',mesh_element(1,cp_en) -!$OMP END CRITICAL (write2out) - call IO_error(650) - return - endif - CPFEM_results(1:3,grain,CPFEM_in,cp_en) = math_RtoEuler(transpose(R))*inDeg ! orientation - CPFEM_results(4 ,grain,CPFEM_in,cp_en) = volfrac ! volume fraction of orientation - enddo ! grain -! - return -! - END SUBROUTINE -! - END MODULE -!############################################################## - + call IO_error(650) + return + endif + CPFEM_results(1:3,g,i,e) = math_RtoEuler(transpose(R))*inDeg ! orientation + CPFEM_results(4 ,g,i,e) = volfrac ! volume fraction of orientation + end do + end do + end do +!$OMP END PARALLEL DO +! + return +! + END SUBROUTINE +! +! +!******************************************************************** +! Initialize crystallite +!******************************************************************** + subroutine crystallite_init() + use mesh, only: mesh_maxNips,mesh_NcpElems + use constitutive, only: constitutive_maxNgrains + + implicit none + + allocate(crystallite_converged(constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)); crystallite_converged = .false. +! +! *** Output to MARC output file *** +!$OMP CRITICAL (write2out) + write(6,*) + write(6,*) 'crystallite Initialization' + write(6,*) + write(6,*) 'crystallite_converged: ', shape(crystallite_converged) + write(6,*) + call flush(6) +!$OMP END CRITICAL (write2out) + return +! + end subroutine +! +! +!******************************************************************** +! Calculates the stress and jacobi (if wanted) for all or a single component +!******************************************************************** + subroutine SingleCrystallite(& + updateJaco,& ! update of Jacobian required + dt,& ! time increment + el_start,& ! first element in element loop + el_end,& ! last element in element loop + CPFEM_in) ! IP number +! + use prec, only: pReal,pInt,pert_Fg,subStepMin, nCutback + use debug + use constitutive + use mesh, only: mesh_element, FE_Nips + use math + use IO, only: IO_error +! use CPFEM + + implicit none +! + logical updateJaco, JacoOK + real(preal) dt + real(pReal), dimension(3,3) :: Fg_pert,Lp_pert, P_pert, Fp_pert, Fe_pert + real(pReal), dimension(6) :: Tstar_v + real(pReal), dimension(constitutive_maxNstatevars) :: state_pert + integer(pInt) el_start, el_end, CPFEM_in, ip_start, ip_end, g, i, e, k, l, iOuter +! + crystallite_converged=.true. +!$OMP PARALLEL DO + do e=el_start,el_end + if(CPFEM_in /= 0_pInt) then + ip_start = CPFEM_in + ip_end = CPFEM_in + else + ip_start = 1 + ip_end = FE_Nips(mesh_element(2,e)) ! iterate over all IPs of this element's type + endif + do i=ip_start,ip_end + do g=1,texture_Ngrains(mesh_element(4,e)) + crystallite_converged(g,i,e)=.false. + end do + end do + end do +!$OMP END PARALLEL DO + constitutive_state_new=constitutive_state_old + CPFEM_Lp_new = CPFEM_Lp_old + iOuter = 0_pInt + do while(any(crystallite_converged(:,:,el_start:el_end))==.false.) +!$OMP PARALLEL DO + do e=el_start,el_end + if(CPFEM_in /= 0_pInt) then + ip_start = CPFEM_in + ip_end = CPFEM_in + else + ip_start = 1 + ip_end = FE_Nips(mesh_element(2,e)) ! iterate over all IPs of this element's type + endif + do i=ip_start,ip_end + do g=1,texture_Ngrains(mesh_element(4,e)) + if(.not.crystallite_converged(g,i,e))& + call IntegrateStress(CPFEM_Tstar_v(:,g,i,e), CPFEM_PK1(:,:,g,i,e), CPFEM_ffn1(:,:,g,i,e),& + CPFEM_Fp_new(:,:,g,i,e), CPFEM_Fe1(:,:,g,i,e), CPFEM_Lp_new(:,:,g,i,e),& + constitutive_state_new(:,g,i,e), dt, g, i, e) + end do + end do + end do +!$OMP END PARALLEL DO +!$OMP PARALLEL DO + do e=el_start,el_end + if(CPFEM_in /= 0_pInt) then + ip_start = CPFEM_in + ip_end = CPFEM_in + else + ip_start = 1 + ip_end = FE_Nips(mesh_element(2,e)) ! iterate over all IPs of this element's type + endif + do i=ip_start,ip_end + do g=1,texture_Ngrains(mesh_element(4,e)) + if(.not.crystallite_converged(g,i,e))& + call UpdateState(CPFEM_Tstar_v(:,g,i,e),constitutive_state_new(:,g,i,e),dt,g,i,e) + end do + end do + end do +!$OMP END PARALLEL DO + iOuter = iOuter + 1_pInt + if (iOuter==Nouter) then +!$OMP CRITICAL (write2out) + write (6,*) 'Terminated outer loop at el,ip,grain',e,i,g +!$OMP CRITICAL (out) + debug_OuterLoopDistribution(iOuter) = debug_OuterLoopDistribution(iOuter)+1 +!$OMP END CRITICAL (out) + call IO_error(600) +!$OMP END CRITICAL (write2out) + endif + end do +!$OMP CRITICAL (out) + debug_OuterLoopDistribution(iOuter) = debug_OuterLoopDistribution(iOuter)+1 +!$OMP END CRITICAL (out) + + if (wantsConstitutiveResults) then ! get the post_results upon request +!$OMP PARALLEL DO + do e=el_start,el_end + if(CPFEM_in /= 0_pInt) then + ip_start = CPFEM_in + ip_end = CPFEM_in + else + ip_start = 1 + ip_end = FE_Nips(mesh_element(2,e)) ! iterate over all IPs of this element's type + endif + do i=ip_start,ip_end + do g=1,texture_Ngrains(mesh_element(4,e)) + CPFEM_results(CPFEM_Nresults+1:CPFEM_Nresults+constitutive_Nresults(g,i,e),g,i,e) =& + constitutive_post_results(CPFEM_Tstar_v(:,g,i,e),constitutive_state_new(:,g,i,e),& + CPFEM_Temperature(i,e),dt,g,i,e) + end do + end do + end do +!$OMP END PARALLEL DO + endif +! +!***** Calculate Jacobian ***** + if(updateJaco) then + if (debugger) then +!$OMP CRITICAL (write2out) + write (6,*) 'Jacobian calc' +!$OMP END CRITICAL (write2out) + endif +! crystallite_converged=.false. +!$OMP PARALLEL DO + do e=el_start,el_end + if(CPFEM_in /= 0_pInt) then + ip_start = CPFEM_in + ip_end = CPFEM_in + else + ip_start = 1 + ip_end = FE_Nips(mesh_element(2,e)) ! iterate over all IPs of this element's type + endif + do i=ip_start,ip_end + do g=1,texture_Ngrains(mesh_element(4,e)) + do k=1,3 + do l=1,3 + crystallite_converged(g,i,e)=.false. + JacoOK=.true. + Fg_pert = CPFEM_ffn1(:,:,g,i,e) ! initialize perturbed Fg + Fg_pert(k,l) = Fg_pert(k,l) + pert_Fg ! perturb single component + Lp_pert = CPFEM_Lp_new(:,:,g,i,e) ! initialize Lp + Fp_pert = CPFEM_Fp_new(:,:,g,i,e) ! initialize Fp + state_pert = constitutive_state_new(:,g,i,e) ! initial guess from end of time step + iOuter=0_pInt + do while(.not.crystallite_converged(g,i,e)) + call IntegrateStress(Tstar_v, P_pert, Fg_pert, Fp_pert, Fe_pert, Lp_pert, state_pert, dt, g, i, e) + call UpdateState(Tstar_v,state_pert,dt,g,i,e) + iOuter = iOuter + 1_pInt + if (iOuter==Nouter) then + JacoOK=.false. + exit + endif + end do +!$OMP CRITICAL (out) + debug_OuterLoopDistribution(iOuter) = debug_OuterLoopDistribution(iOuter)+1 +!$OMP END CRITICAL (out) + if (JacoOK) & + CPFEM_dPdF(:,:,k,l,g,i,e) = (P_pert-CPFEM_PK1(:,:,g,i,e))/pert_Fg ! constructing tangent dP_ij/dFg_kl only if valid forward difference + ! otherwise leave component unchanged + end do + end do + end do + end do + end do +!$OMP END PARALLEL DO + endif +! + return +! + end subroutine +! +!******************************************************************** +! Update the state for a single component +!******************************************************************** + subroutine UpdateState(& + Tstar_v,& ! stress + state,& ! state + dt,& ! time increment + g,& ! grain number + i,& ! integration point number + e& ! element number + ) + use prec, only: pReal,pInt,reltol_Outer + use constitutive, only: constitutive_dotState, constitutive_state_old, constitutive_Nstatevars +! use CPFEM, only: CPFEM_Temperature +! + integer(pInt) g, i, e + real(pReal), dimension(6) :: Tstar_v + real(pReal), dimension(constitutive_Nstatevars(g, i, e)) :: state, ROuter + real(pReal) dt +! + ROuter = state - constitutive_state_old(:,g,i,e) - & + dt*constitutive_dotState(Tstar_v,state,CPFEM_Temperature(i,e),& + g,i,e) ! residuum from evolution of microstructure + state = state - ROuter ! update of microstructure + if (maxval(abs(ROuter/state),state /= 0.0_pReal) < reltol_Outer) crystallite_converged(g,i,e) = .true. +! + return +! + end subroutine +! +! +!******************************************************************** +! Calculates the stress for a single component +!******************************************************************** +!*********************************************************************** +!*** calculation of stress (P), stiffness (dPdF), *** +!*** and announcment of any *** +!*** acceleration of the Newton-Raphson correction *** +!*********************************************************************** + subroutine IntegrateStress(& + Tstar_v,& ! Stress vector + P,& ! first PK stress + Fg_new,& ! new global deformation gradient + Fp_new,& ! new plastic deformation gradient + Fe_new,& ! new "elastic" deformation gradient + Lp,& ! plastic velocity gradient + state_new,& ! new state variable array + dt,& ! time increment + g,& ! grain number + i,& ! integration point number + e) ! element number + +! post_results,& ! plot results from constitutive model +! Fp_new,& ! new plastic deformation gradient +! updateJaco,& ! update of Jacobian required +! Temperature,& ! temperature of crystallite +! Fg_old,& ! old global deformation gradient +! Fp_old,& ! old plastic deformation gradient +! state_old) ! old state variable array +! + use prec, only: pReal,pInt,pert_Fg,subStepMin, nCutback + use debug + use constitutive, only: constitutive_Nstatevars,constitutive_Nresults,constitutive_state_old + use math +! use CPFEM +! + implicit none +! + character(len=128) msg + logical error,success + integer(pInt) e,i,g, nCutbacks, maxCutbacks + real(pReal) Temperature + real(pReal) dt,dt_aim,subFrac,subStep,det + real(pReal), dimension(3,3) :: Lp,Lp_interpolated,inv + real(pReal), dimension(3,3) :: Fg_current,Fg_new,Fg_aim,deltaFg + real(pReal), dimension(3,3) :: Fp_current,Fp_new + real(pReal), dimension(3,3) :: Fe_current,Fe_new + real(pReal), dimension(3,3) :: P + real(pReal), dimension(6) :: Tstar_v + real(pReal), dimension(constitutive_Nstatevars(g,i,e)) :: state_new +! real(pReal), dimension(constitutive_Nstatevars(g,i,e)) :: state_current +! +! debugger= e==1.and.i==1 + deltaFg = Fg_new - CPFEM_ffn(:,:,g,i,e) + subFrac = 0.0_pReal + subStep = 1.0_pReal + nCutbacks = 0_pInt + maxCutbacks = 0_pInt + Fg_current = CPFEM_ffn(:,:,g,i,e) ! initialize to start of inc + Fp_current = CPFEM_Fp_old(:,:,g,i,e) + call math_invert3x3(Fp_current,inv,det,error) + Fe_current = math_mul33x33(Fg_current,inv) +! state_current = state_new + success = .false. ! pretend cutback + dt_aim = 0.0_pReal ! prevent initial Lp interpolation + Temperature=CPFEM_Temperature(i,e) +! +! begin the cutback loop + do while (subStep > subStepMin) ! continue until finished or too much cut backing + if (success) then ! wind forward + Fg_current = Fg_aim + Fe_current = Fe_new + Fp_current = Fp_new +! state_current = state_new + elseif (dt_aim > 0.0_pReal) then + call math_invert3x3(Fg_aim,inv,det,error) ! inv of Fg_aim + Lp_interpolated = 0.5_pReal*Lp + & + 0.5_pReal*(math_I3 - math_mul33x33(Fp_current,& + math_mul33x33(inv,Fe_current)))/dt_aim ! interpolate Lp and L + if (debugger) then +!$OMP CRITICAL (write2out) + write (6,*) 'Lp interpolation' + write (6,'(a,/,3(3(f12.7,x)/))') 'from',Lp(1:3,:) + write (6,'(a,/,3(3(f12.7,x)/))') 'to',Lp_interpolated(1:3,:) +!$OMP END CRITICAL (write2out) + endif + Lp = Lp_interpolated + endif +! + Fg_aim = Fg_current + subStep*deltaFg ! aim for Fg + dt_aim = subStep*dt ! aim for dt + if (debugger) then +!$OMP CRITICAL (write2out) + write (6,*) 'using these values' +! write (6,'(a,/,3(4(f9.3,x)/))') 'state current / MPa',state_current/1e6_pReal + write (6,'(a,/,3(4(f9.3,x)/))') 'state new / MPa',state_new/1e6_pReal + write (6,'(a,/,3(3(f12.7,x)/))') 'Fe current',Fe_current(1:3,:) + write (6,'(a,/,3(3(f12.7,x)/))') 'Fp current',Fp_current(1:3,:) + write (6,'(a,/,3(3(f12.7,x)/))') 'Lp (old=new guess)',Lp(1:3,:) + write (6,'(a20,f,x,a2,x,f)') 'integrating from ',subFrac,'to',(subFrac+subStep) +!$OMP END CRITICAL (write2out) + endif +! + call TimeIntegration(msg,Lp,Fp_new,Fe_new,Tstar_v,P,state_new,dt_aim,e,i,g,Temperature,Fg_aim,Fp_current) +! + + if (msg == 'ok') then + subFrac = subFrac + subStep + subStep = min(1.0_pReal-subFrac, subStep*2.0_pReal) ! accelerate + nCutbacks = 0_pInt ! reset cutback counter + success = .true. ! keep current Lp + else + nCutbacks = nCutbacks + 1 ! record additional cutback + maxCutbacks = max(nCutbacks,maxCutbacks)! remember maximum number of cutbacks + subStep = subStep / 2.0_pReal ! cut time step in half + success = .false. ! force Lp interpolation +! if (debugger) then +!$OMP CRITICAL (write2out) + write (6,*) '>>>>>>>>>>>>>>>>>>>> cutback <<<<<<<<<<<<<<<<<<<<<<' + write (6,*) 'Element, Ip:', e, i + write (6,*) msg +!$OMP END CRITICAL (write2out) +! endif +! + endif + enddo ! potential substepping +! +!$OMP CRITICAL (cutback) + debug_cutbackDistribution(min(nCutback,maxCutbacks)+1) = debug_cutbackDistribution(min(nCutback,maxCutbacks)+1)+1 +!$OMP END CRITICAL (cutback) +! +! debugger = .false. + return + end subroutine + +! +!*********************************************************************** +!*** fully-implicit two-level time integration *** +!*** based on a residuum in Lp and intermediate *** +!*** acceleration of the Newton-Raphson correction *** +!*********************************************************************** + SUBROUTINE TimeIntegration(& + msg,& ! return message + Lpguess,& ! guess of plastic velocity gradient + Fp_new,& ! new plastic deformation gradient + Fe_new,& ! new "elastic" deformation gradient + Tstar_v,& ! Stress vector + P,& ! 1nd PK stress (taken as initial guess if /= 0) + state,& ! current microstructure at end of time inc (taken as guess if /= 0) + dt,& ! time increment + cp_en,& ! element number + ip,& ! integration point number + grain,& ! grain number + Temperature,& ! temperature + Fg_new,& ! new total def gradient + Fp_old) ! former plastic def gradient +! state_current) ! former microstructure + use prec + use debug + use mesh, only: mesh_element + use constitutive, only: constitutive_Nstatevars,constitutive_Microstructure,& + constitutive_homogenizedC,constitutive_LpAndItsTangent + use math + use IO + implicit none +! + character(len=*) msg + logical failed + integer(pInt) cp_en, ip, grain + integer(pInt) iInner,dummy, i,j,k,l,m,n + real(pReal) dt, Temperature, det, p_hydro, leapfrog,maxleap + real(pReal), dimension(6) :: Tstar_v + real(pReal), dimension(9,9) :: dLp,dTdLp,dRdLp,invdRdLp,eye2 + real(pReal), dimension(6,6) :: C_66 + real(pReal), dimension(3,3) :: Fg_new,Fp_new,invFp_new,Fp_old,invFp_old,Fe_new + real(pReal), dimension(3,3) :: P !,Tstar + real(pReal), dimension(3,3) :: Lp,Lpguess,Lpguess_old,Rinner,Rinner_old,A,B,BT,AB,BTA + real(pReal), dimension(3,3,3,3) :: C + real(pReal), dimension(constitutive_Nstatevars(grain, ip, cp_en)) :: state +! + msg = 'ok' ! error-free so far + eye2 = math_identity2nd(9) + + call math_invert3x3(Fp_old,invFp_old,det,failed) ! inversion of Fp_old + if (failed) then + msg = 'inversion Fp_old' + return + endif + + A = math_mul33x33(transpose(invFp_old), math_mul33x33(transpose(Fg_new),math_mul33x33(Fg_new,invFp_old))) +! +! if (all(state == 0.0_pReal)) state = state_current ! former state guessed, if none specified +! iOuter = 0_pInt ! outer counter +! + if (debugger) then +!$OMP CRITICAL (write2out) + write (6,'(a,/,3(3(f12.7,x)/))') 'Fg to be calculated',Fg_new +!$OMP END CRITICAL (write2out) + endif +! +!Outer: do ! outer iteration: State +! iOuter = iOuter+1 +! if (debugger) then +!!$OMP CRITICAL (write2out) +! write (6,'(a,i3)') '---outer ',iOuter +! write (6,'(a,/,3(4(f9.3,x)/))') 'state old / MPa',state_old/1e6_pReal +! write (6,'(a,/,3(4(f9.3,x)/))') 'state / MPa',state/1e6_pReal +! write (6,'(a,/,3(3(f12.7,x)/))') 'Lpguess',Lpguess(1:3,:) +!!$OMP END CRITICAL (write2out) +! endif +! +! if (iOuter > nOuter) then +! msg = 'limit Outer iteration' +!!$OMP CRITICAL (out) +! debug_OuterLoopDistribution(nOuter) = debug_OuterLoopDistribution(nOuter)+1 +!!$OMP END CRITICAL (out) +! return +! endif + call constitutive_Microstructure(state,Temperature,grain,ip,cp_en) + C_66 = constitutive_HomogenizedC(state, grain, ip, cp_en) + C = math_Mandel66to3333(C_66) ! 4th rank elasticity tensor +! + iInner = 0_pInt + leapfrog = 1.0_pReal ! correction as suggested by invdRdLp-step + maxleap = 1024.0_pReal ! preassign maximum acceleration level +! + Lpguess_old = Lpguess ! consider present Lpguess good +! +Inner: do ! inner iteration: Lp + iInner = iInner+1 +! if (debugger) then +!!$OMP CRITICAL (write2out) +! write (6,'(a,i3)') 'inner ',iInner +! if (iInner < 3) then +! write (6,'(a,/,3(3(f12.7,x)/))') 'Lpguess',Lpguess(1:3,:) +! endif +!!$OMP END CRITICAL (write2out) +! endif + if (iInner > nInner) then ! too many loops required + Lpguess = Lpguess_old ! do not trust the last update but resort to former one + msg = 'limit Inner iteration' +!$OMP CRITICAL (in) + debug_InnerLoopDistribution(nInner) = debug_InnerLoopDistribution(nInner)+1 +!$OMP END CRITICAL (in) + return + endif +! + B = math_i3 - dt*Lpguess + BT = transpose(B) + AB = math_mul33x33(A,B) + BTA = math_mul33x33(BT,A) + Tstar_v = 0.5_pReal*math_mul66x6(C_66,math_mandel33to6(math_mul33x33(BT,AB)-math_I3)) +! Tstar = math_Mandel6to33(Tstar_v) + p_hydro=(Tstar_v(1)+Tstar_v(2)+Tstar_v(3))/3.0_pReal + forall(i=1:3) Tstar_v(i) = Tstar_v(i)-p_hydro ! subtract hydrostatic pressure + call constitutive_LpAndItsTangent(Lp,dLp, & + Tstar_v,state,Temperature,grain,ip,cp_en) +! + Rinner = Lpguess - Lp ! update current residuum +! + if (.not.(any(Rinner/=Rinner)) .and. & ! exclude any NaN in residuum + ( (maxval(abs(Rinner)) < abstol_Inner) .or. & ! below abs tol .or. + ( any(abs(dt*Lpguess) > relevantStrain) .and. & ! worth checking? .and. + maxval(abs(Rinner/Lpguess),abs(dt*Lpguess) > relevantStrain) < reltol_Inner & ! below rel tol + ) & + ) & + ) & + exit Inner ! convergence +! +! check for acceleration/deceleration in Newton--Raphson correction +! + if (any(Rinner/=Rinner) .and. & ! NaN occured at regular speed + leapfrog == 1.0) then + Lpguess = Lpguess_old ! restore known good guess + msg = 'NaN present' ! croak for cutback + return + + elseif (leapfrog > 1.0_pReal .and. & ! at fast pace ? + (sum(Rinner*Rinner) > sum(Rinner_old*Rinner_old) .or. & ! worse residuum + sum(Rinner*Rinner_old) < 0.0_pReal) .or. & ! residuum changed sign (overshoot) + any(Rinner/=Rinner) ) then ! NaN + maxleap = 0.5_pReal * leapfrog ! limit next acceleration + leapfrog = 1.0_pReal ! grinding halt + + else ! better residuum + dTdLp = 0.0_pReal ! calc dT/dLp + forall (i=1:3,j=1:3,k=1:3,l=1:3,m=1:3,n=1:3) & + dTdLp(3*(i-1)+j,3*(k-1)+l) = dTdLp(3*(i-1)+j,3*(k-1)+l) + & + C(i,j,l,n)*AB(k,n)+C(i,j,m,l)*BTA(m,k) + dTdLp = -0.5_pReal*dt*dTdLp + dRdLp = eye2 - math_mul99x99(dLp,dTdLp) ! calc dR/dLp + invdRdLp = 0.0_pReal + call math_invert(9,dRdLp,invdRdLp,dummy,failed) ! invert dR/dLp --> dLp/dR + if (failed) then + msg = 'inversion dR/dLp' + if (debugger) then +!$OMP CRITICAL (write2out) + write (6,*) msg + write (6,'(a,/,9(9(e9.3,x)/))') 'dRdLp', dRdLp(1:9,:) + write (6,'(a,/,3(4(f9.3,x)/))') 'state / MPa',state/1e6_pReal + write (6,'(a,/,3(3(f12.7,x)/))') 'Lpguess',Lpguess(1:3,:) + write (6,'(a,/,3(3(e12.7,x)/))') 'Lp',Lp(1:3,:) + write (6,'(a,/,6(f9.3,x))') 'Tstar / MPa',Tstar_v/1e6_pReal +!$OMP END CRITICAL (write2out) + endif + return + endif +! + Rinner_old = Rinner ! remember current residuum + Lpguess_old = Lpguess ! remember current Lp guess + if (iInner > 1 .and. leapfrog < maxleap) & + leapfrog = 2.0_pReal * leapfrog ! accelerate if ok + endif +! + Lpguess = Lpguess_old ! start from current guess + Rinner = Rinner_old ! use current residuum + forall (i=1:3,j=1:3,k=1:3,l=1:3) & ! leapfrog to updated Lpguess + Lpguess(i,j) = Lpguess(i,j) - leapfrog*invdRdLp(3*(i-1)+j,3*(k-1)+l)*Rinner(k,l) + enddo Inner +! +!$OMP CRITICAL (in) + debug_InnerLoopDistribution(iInner) = debug_InnerLoopDistribution(iInner)+1 +!$OMP END CRITICAL (in) +! ROuter = state - state_old - & +! dt*constitutive_dotState(Tstar_v,state,Temperature,& +! grain,ip,cp_en) ! residuum from evolution of microstructure +! state = state - ROuter ! update of microstructure +! +! if (iOuter==nOuter) then +!!$OMP CRITICAL (write2out) +! write (6,*) 'Terminated outer loop at el,ip,grain',cp_en,ip,grain +!!$OMP END CRITICAL (write2out) +! exit Outer +! endif +! if (maxval(abs(Router/state),state /= 0.0_pReal) < reltol_Outer) exit Outer +! enddo Outer +! +!!$OMP CRITICAL (out) +! debug_OuterLoopDistribution(iOuter) = debug_OuterLoopDistribution(iOuter)+1 +!!$OMP END CRITICAL (out) + invFp_new = math_mul33x33(invFp_old,B) + call math_invert3x3(invFp_new,Fp_new,det,failed) + if (failed) then + msg = 'inversion Fp_new^-1' + return + endif +! +! if (wantsConstitutiveResults) then ! get the post_results upon request +! results = 0.0_pReal +! results = constitutive_post_results(Tstar_v,state,Temperature,dt,grain,ip,cp_en) +! endif +! + Fp_new = Fp_new*det**(1.0_pReal/3.0_pReal) ! regularize Fp by det = det(InvFp_new) !! + forall (i=1:3) Tstar_v(i) = Tstar_v(i)+p_hydro ! add hydrostatic component back + Fe_new = math_mul33x33(Fg_new,invFp_new) ! calc resulting Fe +! P = math_mul33x33(Fe_new,math_mul33x33(Tstar,transpose(invFp_new))) ! first PK stress + P = math_mul33x33(Fe_new,math_mul33x33(math_Mandel6to33(Tstar_v),transpose(invFp_new))) ! first PK stress + + return +! + END SUBROUTINE +! + END MODULE +!############################################################## + diff --git a/trunk/CPFEM_Taylor_sequential.f90 b/trunk/CPFEM_Taylor_sequential.f90 index 41d523bf3..11e1e22b7 100644 --- a/trunk/CPFEM_Taylor_sequential.f90 +++ b/trunk/CPFEM_Taylor_sequential.f90 @@ -9,23 +9,34 @@ ! **************************************************************** ! *** General variables for the material behaviour calculation *** ! **************************************************************** - real(pReal), dimension (:,:), allocatable :: CPFEM_Temperature - real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_ffn_bar - real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_ffn1_bar - real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_PK1_bar - real(pReal), dimension (:,:,:,:,:,:),allocatable :: CPFEM_dPdF_bar - 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 - real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_Fp_old - real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_Fp_new + 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_Fe1 + real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_Tstar_v real(pReal), parameter :: CPFEM_odd_stress = 1e15_pReal, CPFEM_odd_jacobian = 1e50_pReal integer(pInt) :: CPFEM_Nresults = 4_pInt ! three Euler angles plus volume fraction logical :: CPFEM_init_done = .false. ! remember whether init has been done already logical :: CPFEM_calc_done = .false. ! remember whether first IP has already calced the results logical :: CPFEM_results_aged = .false. ! remember whether results have been aged at inc start +! *** Solution at single crystallite level *** +! + logical, dimension (:,:,:),allocatable :: crystallite_converged !individual covergence flag per grain ! CONTAINS ! @@ -46,12 +57,18 @@ 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)) + 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_ffn1_bar (3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_ffn1_bar = CPFEM_ffn_bar - allocate(CPFEM_PK1_bar (3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_PK1_bar = 0.0_pReal + allocate(CPFEM_ffn(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) + forall(g=1:constitutive_maxNgrains,e=1:mesh_NcpElems,i=1:mesh_maxNips) CPFEM_ffn(:,:,g,i,e) = math_I3 + allocate(CPFEM_ffn1_bar(3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_ffn1_bar = CPFEM_ffn_bar + allocate(CPFEM_ffn1(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_ffn1 = CPFEM_ffn + allocate(CPFEM_PK1_bar(3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_PK1_bar = 0.0_pReal + allocate(CPFEM_PK1(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_PK1 = 0.0_pReal allocate(CPFEM_dPdF_bar(3,3,3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dPdF_bar = 0.0_pReal + allocate(CPFEM_dPdF_bar_old(3,3,3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dPdF_bar_old = 0.0_pReal + allocate(CPFEM_dPdF(3,3,3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dPdF = 0.0_pReal allocate(CPFEM_stress_bar(6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_stress_bar = 0.0_pReal allocate(CPFEM_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 @@ -61,31 +78,44 @@ CPFEM_results = 0.0_pReal ! ! *** Plastic velocity gradient *** - allocate(CPFEM_Lp(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Lp = 0.0_pReal + allocate(CPFEM_Lp_old(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Lp_old = 0.0_pReal + allocate(CPFEM_Lp_new(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Lp_new = 0.0_pReal ! *** Plastic deformation gradient at (t=t0) and (t=t1) *** allocate(CPFEM_Fp_new(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Fp_new = 0.0_pReal allocate(CPFEM_Fp_old(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) forall (e=1:mesh_NcpElems,i=1:mesh_maxNips,g=1:constitutive_maxNgrains) & CPFEM_Fp_old(:,:,g,i,e) = math_EulerToR(constitutive_EulerAngles(:,g,i,e)) ! plastic def gradient reflects init orientation +! *** Elastic deformation gradient at (t=t1) *** + allocate(CPFEM_Fe1(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Fe1 = 0.0_pReal +! *** Stress vector at (t=t1) *** + allocate(CPFEM_Tstar_v(6,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Tstar_v = 0.0_pReal ! ! *** Output to MARC output file *** !$OMP CRITICAL (write2out) write(6,*) write(6,*) 'CPFEM Initialization' write(6,*) - write(6,*) 'CPFEM_Temperature: ', shape(CPFEM_Temperature) - write(6,*) 'CPFEM_ffn_bar: ', shape(CPFEM_ffn_bar) - write(6,*) 'CPFEM_ffn1_bar: ', shape(CPFEM_ffn1_bar) - write(6,*) 'CPFEM_PK1_bar: ', shape(CPFEM_PK1_bar) - write(6,*) 'CPFEM_dPdF_bar: ', shape(CPFEM_dPdF_bar) - write(6,*) 'CPFEM_stress_bar: ', shape(CPFEM_stress_bar) - write(6,*) 'CPFEM_jaco_bar: ', shape(CPFEM_jaco_bar) + 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: ', shape(CPFEM_Lp) - write(6,*) 'CPFEM_Fp_old: ', shape(CPFEM_Fp_old) - write(6,*) 'CPFEM_Fp_new: ', shape(CPFEM_Fp_new) + 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_Fe1: ', shape(CPFEM_Fe1) + write(6,*) 'CPFEM_Tstar_v: ', shape(CPFEM_Tstar_v) write(6,*) call flush(6) !$OMP END CRITICAL (write2out) @@ -142,6 +172,7 @@ call mesh_init() call lattice_init() call constitutive_init() + call crystallite_init() call CPFEM_init(Temperature) CPFEM_init_done = .true. endif @@ -186,9 +217,9 @@ math_I3(j,l)*CPFEM_ffn1_bar(i,m,CPFEM_in,cp_en)*CPFEM_PK1_bar(k,m,CPFEM_in,cp_en) + & 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)) - 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) + 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) case (3) ! *** collect and return odd result *** CPFEM_Temperature(CPFEM_in,cp_en) = Temperature @@ -217,111 +248,704 @@ END SUBROUTINE ! ! -!********************************************************** -!*** calculate the material point behaviour *** -!********************************************************** - SUBROUTINE CPFEM_MaterialPoint(& - updateJaco,& ! flag to initiate Jacobian updating - CPFEM_dt,& ! Time increment (dt) - CPFEM_in,& ! Integration point number - cp_en) ! Element number -! - use prec - use FEsolving, only: theCycle - use debug - use math, only: math_pDecomposition,math_RtoEuler,inDeg,math_I3,math_invert3x3,math_permut,math_invert,math_delta - use IO, only: IO_error - use mesh, only: mesh_element - use crystallite - use constitutive - implicit none -! - character(len=128) msg - integer(pInt) cp_en,CPFEM_in,grain - logical updateJaco,error - real(pReal) CPFEM_dt,volfrac - real(pReal), dimension(3,3) :: U,R,Fe1 - real(pReal), dimension(3,3) :: PK1 - real(pReal), dimension(3,3,3,3) :: dPdF,dPdF_bar_old -! - CPFEM_PK1_bar(:,:,CPFEM_in,cp_en) = 0.0_pReal ! zero out average first PK stress - if (updateJaco) then - dPdF_bar_old = CPFEM_dPdF_bar(:,:,:,:,CPFEM_in,cp_en) ! remember former average consistent tangent - CPFEM_dPdF_bar(:,:,:,:,CPFEM_in,cp_en) = 0.0_pReal ! zero out avg consistent tangent for later assembly - endif - - do grain = 1,texture_Ngrains(mesh_element(4,cp_en)) - dPdF = dPdF_bar_old ! preguess consistent tangent of grain with avg - - if (debugger) then -!$OMP CRITICAL (write2out) - write (6,*) 'single crystallite integrating.',cp_en,CPFEM_in,grain - write (6,'(a,/,3(3(f12.7,x)/))') 'Fg',CPFEM_ffn1_bar(:,:,CPFEM_in,cp_en) - write (6,'(a,/,3(3(f12.7,x)/))') 'Lp (guess)',CPFEM_Lp(1:3,:,grain,CPFEM_in,cp_en) - write (6,'(a,/,3(3(f12.7,x)/))') 'Fp (old)',CPFEM_Fp_old(1:3,:,grain,CPFEM_in,cp_en) - write (6,'(a,/,3(4(f9.3,x)/))') 'state (old) / MPa',constitutive_state_old(:,grain,CPFEM_in,cp_en)/1e6_pReal - write (6,'(a,/,3(4(f9.3,x)/))') 'state (new) / MPa',constitutive_state_new(:,grain,CPFEM_in,cp_en)/1e6_pReal - write (6,*) -!$OMP END CRITICAL (write2out) - endif - - call SingleCrystallite(msg,PK1,dPdF,& - CPFEM_results(CPFEM_Nresults+1:CPFEM_Nresults+constitutive_Nresults(grain,CPFEM_in,cp_en),& - grain,CPFEM_in,cp_en),& - CPFEM_Lp(:,:,grain,CPFEM_in,cp_en),& - CPFEM_Fp_new(:,:,grain,CPFEM_in,cp_en),Fe1,constitutive_state_new(:,grain,CPFEM_in,cp_en),& ! output up to here - CPFEM_dt,cp_en,CPFEM_in,grain,updateJaco,& - CPFEM_Temperature(CPFEM_in,cp_en),& - CPFEM_ffn1_bar(:,:,CPFEM_in,cp_en),CPFEM_ffn_bar(:,:,CPFEM_in,cp_en),& - CPFEM_Fp_old(:,:,grain,CPFEM_in,cp_en),constitutive_state_old(:,grain,CPFEM_in,cp_en)) - - if (msg /= 'ok') then ! solution not reached --> exit -!$OMP CRITICAL (write2out) - write(6,*) 'grain loop failed to converge @ EL:',cp_en,' IP:',CPFEM_in -!$OMP END CRITICAL (write2out) - call IO_error(600) - return - endif - - if (debugger) then -!$OMP CRITICAL (write2out) - write (6,*) msg - write (6,*) 'single crystallite convergence reached.',cp_en,CPFEM_in,grain - write (6,'(a,/,3(3(f12.7,x)/))') 'Lp',CPFEM_Lp(1:3,:,grain,CPFEM_in,cp_en) - write (6,'(a,/,3(3(f12.7,x)/))') 'Fp (new)',CPFEM_Fp_new(1:3,:,grain,CPFEM_in,cp_en) - write (6,'(a,/,3(4(f9.3,x)/))') 'state (new)/ MPa',constitutive_state_new(:,grain,CPFEM_in,cp_en)/1e6_pReal - write (6,'(a,/,3(3(f9.3,x)/))') 'P / MPa',PK1/1e6_pReal - write (6,'(a,/,9(9(f9.3,x)/))') 'dP/dF / GPa',dPdF/1e9_pReal -!$OMP END CRITICAL (write2out) - endif - - volfrac = constitutive_matVolFrac(grain,CPFEM_in,cp_en)*constitutive_texVolFrac(grain,CPFEM_in,cp_en) - CPFEM_PK1_bar(:,:,CPFEM_in,cp_en) = CPFEM_PK1_bar(:,:,CPFEM_in,cp_en) + volfrac*PK1 - if (updateJaco) CPFEM_dPdF_bar(:,:,:,:,CPFEM_in,cp_en) = & - CPFEM_dPdF_bar(:,:,:,:,CPFEM_in,cp_en) + volfrac*dPdF ! add up crystallite stiffnesses - ! (may have "holes" corresponding - ! to former avg tangent) -! -! update results plotted in MENTAT - call math_pDecomposition(Fe1,U,R,error) ! polar decomposition - if (error) then -!$OMP CRITICAL (write2out) - write(6,*) 'polar decomposition of', Fe1 - write(6,*) 'Grain: ',grain - write(6,*) 'Integration point: ',CPFEM_in - write(6,*) 'Element: ',mesh_element(1,cp_en) -!$OMP END CRITICAL (write2out) - call IO_error(650) - return - endif - CPFEM_results(1:3,grain,CPFEM_in,cp_en) = math_RtoEuler(transpose(R))*inDeg ! orientation - CPFEM_results(4 ,grain,CPFEM_in,cp_en) = volfrac ! volume fraction of orientation - enddo ! grain -! - return -! - END SUBROUTINE -! - END MODULE -!############################################################## - +!********************************************************** +!*** calculate the material point behaviour *** +!********************************************************** + SUBROUTINE CPFEM_MaterialPoint(& + updateJaco,& ! flag to initiate Jacobian updating + CPFEM_dt,& ! Time increment (dt) + CPFEM_in,& ! Integration point number + cp_en) ! Element number +! + use prec + use FEsolving, only: theCycle + use debug + use math, only: math_pDecomposition,math_RtoEuler,inDeg,math_I3,math_invert3x3,math_permut,math_invert,math_delta + use IO, only: IO_error + use mesh, only: mesh_element, mesh_NcpElems, FE_Nips +! use crystallite + use constitutive + implicit none +! + integer(pInt) cp_en,CPFEM_in,g,i,e + integer(pInt) el_start, el_end, ip_start, ip_end + logical updateJaco,error + real(pReal) CPFEM_dt,volfrac + real(pReal), dimension(3,3) :: U,R !,Fe1 +! real(pReal), dimension(3,3) :: PK1 +! real(pReal), dimension(3,3,3,3) :: dPdF,dPdF_bar_old +! + CPFEM_PK1_bar = 0.0_pReal ! zero out average first PK stress +!initialize element loop + if (cp_en /= 0_pInt) then + el_start = cp_en + el_end = cp_en + else + el_start = 1_pInt + el_end = mesh_NcpElems + endif +! prescribe FFN and FFN1 depending on homogenization scheme +!$OMP PARALLEL DO + do e=el_start,el_end + if(CPFEM_in /= 0_pInt) then + ip_start = CPFEM_in + ip_end = CPFEM_in + else + ip_start = 1 + ip_end = FE_Nips(mesh_element(2,e)) ! iterate over all IPs of this element's type + endif + do i=ip_start,ip_end + do g=1,texture_Ngrains(mesh_element(4,e)) + CPFEM_ffn(:,:,g,i,e) = CPFEM_ffn_bar(:,:,i,e) !Taylor homogenization + CPFEM_ffn1(:,:,g,i,e) = CPFEM_ffn1_bar(:,:,i,e) !Taylor homogenization + end do + end do + end do +!$OMP END PARALLEL DO +! calculate stress, update state and update jacobian in case needed for all or one ip + if (updateJaco) then + CPFEM_dPdF_bar_old = CPFEM_dPdF_bar ! remember former average consistent tangent + CPFEM_dPdF_bar = 0.0_pReal ! zero out avg consistent tangent for later assembly + endif + call SingleCrystallite(updateJaco,CPFEM_dt,el_start,el_end,CPFEM_in) +!****************************************************************************************************** +! check convergence of homogenization in case needed +!****************************************************************************************************** +! calculate average quantities per ip and post results +!$OMP PARALLEL DO + do e=el_start,el_end + if(CPFEM_in /= 0_pInt) then + ip_start = CPFEM_in + ip_end = CPFEM_in + else + ip_start = 1 + ip_end = FE_Nips(mesh_element(2,e)) ! iterate over all IPs of this element's type + endif + do i=ip_start,ip_end + do g=1,texture_Ngrains(mesh_element(4,e)) + volfrac = constitutive_matVolFrac(g,i,e)*constitutive_texVolFrac(g,i,e) + CPFEM_PK1_bar(:,:,i,e) = CPFEM_PK1_bar(:,:,i,e) + volfrac * CPFEM_PK1(:,:,g,i,e) + if (updateJaco) CPFEM_dPdF_bar(:,:,:,:,i,e) = & + CPFEM_dPdF_bar(:,:,:,:,i,e) + volfrac * CPFEM_dPdF(:,:,:,:,g,i,e) ! add up crystallite stiffnesses + ! (may have "holes" corresponding + ! to former avg tangent) +! update results plotted in MENTAT + call math_pDecomposition(CPFEM_Fe1(:,:,g,i,e),U,R,error) ! polar decomposition + if (error) then +!$OMP CRITICAL (write2out) + write(6,*) 'polar decomposition of', CPFEM_Fe1(:,:,g,i,e) + write(6,*) 'Grain: ',g + write(6,*) 'Integration point: ',i + write(6,*) 'Element: ',mesh_element(1,e) +!$OMP END CRITICAL (write2out) + call IO_error(650) + return + endif + CPFEM_results(1:3,g,i,e) = math_RtoEuler(transpose(R))*inDeg ! orientation + CPFEM_results(4 ,g,i,e) = volfrac ! volume fraction of orientation + end do + end do + end do +!$OMP END PARALLEL DO +! + return +! + END SUBROUTINE +! +! +!******************************************************************** +! Initialize crystallite +!******************************************************************** + subroutine crystallite_init() + use mesh, only: mesh_maxNips,mesh_NcpElems + use constitutive, only: constitutive_maxNgrains + + implicit none + + allocate(crystallite_converged(constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)); crystallite_converged = .false. +! +! *** Output to MARC output file *** +!$OMP CRITICAL (write2out) + write(6,*) + write(6,*) 'crystallite Initialization' + write(6,*) + write(6,*) 'crystallite_converged: ', shape(crystallite_converged) + write(6,*) + call flush(6) +!$OMP END CRITICAL (write2out) + return +! + end subroutine +! +! +!******************************************************************** +! Calculates the stress and jacobi (if wanted) for all or a single component +!******************************************************************** + subroutine SingleCrystallite(& + updateJaco,& ! update of Jacobian required + dt,& ! time increment + el_start,& ! first element in element loop + el_end,& ! last element in element loop + CPFEM_in) ! IP number +! + use prec, only: pReal,pInt,pert_Fg,subStepMin, nCutback + use debug + use constitutive + use mesh, only: mesh_element, FE_Nips + use math + use IO, only: IO_error +! use CPFEM + + implicit none +! + logical updateJaco, JacoOK + real(preal) dt + real(pReal), dimension(3,3) :: Fg_pert,Lp_pert, P_pert, Fp_pert, Fe_pert + real(pReal), dimension(6) :: Tstar_v + real(pReal), dimension(constitutive_maxNstatevars) :: state_pert + integer(pInt) el_start, el_end, CPFEM_in, ip_start, ip_end, g, i, e, k, l, iOuter +! + crystallite_converged=.true. +!$OMP PARALLEL DO + do e=el_start,el_end + if(CPFEM_in /= 0_pInt) then + ip_start = CPFEM_in + ip_end = CPFEM_in + else + ip_start = 1 + ip_end = FE_Nips(mesh_element(2,e)) ! iterate over all IPs of this element's type + endif + do i=ip_start,ip_end + do g=1,texture_Ngrains(mesh_element(4,e)) + crystallite_converged(g,i,e)=.false. + end do + end do + end do +!$OMP END PARALLEL DO + constitutive_state_new=constitutive_state_old + CPFEM_Lp_new = CPFEM_Lp_old + iOuter = 0_pInt + do while(any(crystallite_converged(:,:,el_start:el_end))==.false.) +!$OMP PARALLEL DO + do e=el_start,el_end + if(CPFEM_in /= 0_pInt) then + ip_start = CPFEM_in + ip_end = CPFEM_in + else + ip_start = 1 + ip_end = FE_Nips(mesh_element(2,e)) ! iterate over all IPs of this element's type + endif + do i=ip_start,ip_end + do g=1,texture_Ngrains(mesh_element(4,e)) + if(.not.crystallite_converged(g,i,e))& + call IntegrateStress(CPFEM_Tstar_v(:,g,i,e), CPFEM_PK1(:,:,g,i,e), CPFEM_ffn1(:,:,g,i,e),& + CPFEM_Fp_new(:,:,g,i,e), CPFEM_Fe1(:,:,g,i,e), CPFEM_Lp_new(:,:,g,i,e),& + constitutive_state_new(:,g,i,e), dt, g, i, e) + end do + end do + end do +!$OMP END PARALLEL DO +!$OMP PARALLEL DO + do e=el_start,el_end + if(CPFEM_in /= 0_pInt) then + ip_start = CPFEM_in + ip_end = CPFEM_in + else + ip_start = 1 + ip_end = FE_Nips(mesh_element(2,e)) ! iterate over all IPs of this element's type + endif + do i=ip_start,ip_end + do g=1,texture_Ngrains(mesh_element(4,e)) + if(.not.crystallite_converged(g,i,e))& + call UpdateState(CPFEM_Tstar_v(:,g,i,e),constitutive_state_new(:,g,i,e),dt,g,i,e) + end do + end do + end do +!$OMP END PARALLEL DO + iOuter = iOuter + 1_pInt + if (iOuter==Nouter) then +!$OMP CRITICAL (write2out) + write (6,*) 'Terminated outer loop at el,ip,grain',e,i,g +!$OMP CRITICAL (out) + debug_OuterLoopDistribution(iOuter) = debug_OuterLoopDistribution(iOuter)+1 +!$OMP END CRITICAL (out) + call IO_error(600) +!$OMP END CRITICAL (write2out) + endif + end do +!$OMP CRITICAL (out) + debug_OuterLoopDistribution(iOuter) = debug_OuterLoopDistribution(iOuter)+1 +!$OMP END CRITICAL (out) + + if (wantsConstitutiveResults) then ! get the post_results upon request +!$OMP PARALLEL DO + do e=el_start,el_end + if(CPFEM_in /= 0_pInt) then + ip_start = CPFEM_in + ip_end = CPFEM_in + else + ip_start = 1 + ip_end = FE_Nips(mesh_element(2,e)) ! iterate over all IPs of this element's type + endif + do i=ip_start,ip_end + do g=1,texture_Ngrains(mesh_element(4,e)) + CPFEM_results(CPFEM_Nresults+1:CPFEM_Nresults+constitutive_Nresults(g,i,e),g,i,e) =& + constitutive_post_results(CPFEM_Tstar_v(:,g,i,e),constitutive_state_new(:,g,i,e),& + CPFEM_Temperature(i,e),dt,g,i,e) + end do + end do + end do +!$OMP END PARALLEL DO + endif +! +!***** Calculate Jacobian ***** + if(updateJaco) then + if (debugger) then +!$OMP CRITICAL (write2out) + write (6,*) 'Jacobian calc' +!$OMP END CRITICAL (write2out) + endif +! crystallite_converged=.false. +!$OMP PARALLEL DO + do e=el_start,el_end + if(CPFEM_in /= 0_pInt) then + ip_start = CPFEM_in + ip_end = CPFEM_in + else + ip_start = 1 + ip_end = FE_Nips(mesh_element(2,e)) ! iterate over all IPs of this element's type + endif + do i=ip_start,ip_end + do g=1,texture_Ngrains(mesh_element(4,e)) + do k=1,3 + do l=1,3 + crystallite_converged(g,i,e)=.false. + JacoOK=.true. + Fg_pert = CPFEM_ffn1(:,:,g,i,e) ! initialize perturbed Fg + Fg_pert(k,l) = Fg_pert(k,l) + pert_Fg ! perturb single component + Lp_pert = CPFEM_Lp_new(:,:,g,i,e) ! initialize Lp + Fp_pert = CPFEM_Fp_new(:,:,g,i,e) ! initialize Fp + state_pert = constitutive_state_new(:,g,i,e) ! initial guess from end of time step + iOuter=0_pInt + do while(.not.crystallite_converged(g,i,e)) + call IntegrateStress(Tstar_v, P_pert, Fg_pert, Fp_pert, Fe_pert, Lp_pert, state_pert, dt, g, i, e) + call UpdateState(Tstar_v,state_pert,dt,g,i,e) + iOuter = iOuter + 1_pInt + if (iOuter==Nouter) then + JacoOK=.false. + exit + endif + end do +!$OMP CRITICAL (out) + debug_OuterLoopDistribution(iOuter) = debug_OuterLoopDistribution(iOuter)+1 +!$OMP END CRITICAL (out) + if (JacoOK) & + CPFEM_dPdF(:,:,k,l,g,i,e) = (P_pert-CPFEM_PK1(:,:,g,i,e))/pert_Fg ! constructing tangent dP_ij/dFg_kl only if valid forward difference + ! otherwise leave component unchanged + end do + end do + end do + end do + end do +!$OMP END PARALLEL DO + endif +! + return +! + end subroutine +! +!******************************************************************** +! Update the state for a single component +!******************************************************************** + subroutine UpdateState(& + Tstar_v,& ! stress + state,& ! state + dt,& ! time increment + g,& ! grain number + i,& ! integration point number + e& ! element number + ) + use prec, only: pReal,pInt,reltol_Outer + use constitutive, only: constitutive_dotState, constitutive_state_old, constitutive_Nstatevars +! use CPFEM, only: CPFEM_Temperature +! + integer(pInt) g, i, e + real(pReal), dimension(6) :: Tstar_v + real(pReal), dimension(constitutive_Nstatevars(g, i, e)) :: state, ROuter + real(pReal) dt +! + ROuter = state - constitutive_state_old(:,g,i,e) - & + dt*constitutive_dotState(Tstar_v,state,CPFEM_Temperature(i,e),& + g,i,e) ! residuum from evolution of microstructure + state = state - ROuter ! update of microstructure + if (maxval(abs(ROuter/state),state /= 0.0_pReal) < reltol_Outer) crystallite_converged(g,i,e) = .true. +! + return +! + end subroutine +! +! +!******************************************************************** +! Calculates the stress for a single component +!******************************************************************** +!*********************************************************************** +!*** calculation of stress (P), stiffness (dPdF), *** +!*** and announcment of any *** +!*** acceleration of the Newton-Raphson correction *** +!*********************************************************************** + subroutine IntegrateStress(& + Tstar_v,& ! Stress vector + P,& ! first PK stress + Fg_new,& ! new global deformation gradient + Fp_new,& ! new plastic deformation gradient + Fe_new,& ! new "elastic" deformation gradient + Lp,& ! plastic velocity gradient + state_new,& ! new state variable array + dt,& ! time increment + g,& ! grain number + i,& ! integration point number + e) ! element number + +! post_results,& ! plot results from constitutive model +! Fp_new,& ! new plastic deformation gradient +! updateJaco,& ! update of Jacobian required +! Temperature,& ! temperature of crystallite +! Fg_old,& ! old global deformation gradient +! Fp_old,& ! old plastic deformation gradient +! state_old) ! old state variable array +! + use prec, only: pReal,pInt,pert_Fg,subStepMin, nCutback + use debug + use constitutive, only: constitutive_Nstatevars,constitutive_Nresults,constitutive_state_old + use math +! use CPFEM +! + implicit none +! + character(len=128) msg + logical error,success + integer(pInt) e,i,g, nCutbacks, maxCutbacks + real(pReal) Temperature + real(pReal) dt,dt_aim,subFrac,subStep,det + real(pReal), dimension(3,3) :: Lp,Lp_interpolated,inv + real(pReal), dimension(3,3) :: Fg_current,Fg_new,Fg_aim,deltaFg + real(pReal), dimension(3,3) :: Fp_current,Fp_new + real(pReal), dimension(3,3) :: Fe_current,Fe_new + real(pReal), dimension(3,3) :: P + real(pReal), dimension(6) :: Tstar_v + real(pReal), dimension(constitutive_Nstatevars(g,i,e)) :: state_new +! real(pReal), dimension(constitutive_Nstatevars(g,i,e)) :: state_current +! +! debugger= e==1.and.i==1 + deltaFg = Fg_new - CPFEM_ffn(:,:,g,i,e) + subFrac = 0.0_pReal + subStep = 1.0_pReal + nCutbacks = 0_pInt + maxCutbacks = 0_pInt + Fg_current = CPFEM_ffn(:,:,g,i,e) ! initialize to start of inc + Fp_current = CPFEM_Fp_old(:,:,g,i,e) + call math_invert3x3(Fp_current,inv,det,error) + Fe_current = math_mul33x33(Fg_current,inv) +! state_current = state_new + success = .false. ! pretend cutback + dt_aim = 0.0_pReal ! prevent initial Lp interpolation + Temperature=CPFEM_Temperature(i,e) +! +! begin the cutback loop + do while (subStep > subStepMin) ! continue until finished or too much cut backing + if (success) then ! wind forward + Fg_current = Fg_aim + Fe_current = Fe_new + Fp_current = Fp_new +! state_current = state_new + elseif (dt_aim > 0.0_pReal) then + call math_invert3x3(Fg_aim,inv,det,error) ! inv of Fg_aim + Lp_interpolated = 0.5_pReal*Lp + & + 0.5_pReal*(math_I3 - math_mul33x33(Fp_current,& + math_mul33x33(inv,Fe_current)))/dt_aim ! interpolate Lp and L + if (debugger) then +!$OMP CRITICAL (write2out) + write (6,*) 'Lp interpolation' + write (6,'(a,/,3(3(f12.7,x)/))') 'from',Lp(1:3,:) + write (6,'(a,/,3(3(f12.7,x)/))') 'to',Lp_interpolated(1:3,:) +!$OMP END CRITICAL (write2out) + endif + Lp = Lp_interpolated + endif +! + Fg_aim = Fg_current + subStep*deltaFg ! aim for Fg + dt_aim = subStep*dt ! aim for dt + if (debugger) then +!$OMP CRITICAL (write2out) + write (6,*) 'using these values' +! write (6,'(a,/,3(4(f9.3,x)/))') 'state current / MPa',state_current/1e6_pReal + write (6,'(a,/,3(4(f9.3,x)/))') 'state new / MPa',state_new/1e6_pReal + write (6,'(a,/,3(3(f12.7,x)/))') 'Fe current',Fe_current(1:3,:) + write (6,'(a,/,3(3(f12.7,x)/))') 'Fp current',Fp_current(1:3,:) + write (6,'(a,/,3(3(f12.7,x)/))') 'Lp (old=new guess)',Lp(1:3,:) + write (6,'(a20,f,x,a2,x,f)') 'integrating from ',subFrac,'to',(subFrac+subStep) +!$OMP END CRITICAL (write2out) + endif +! + call TimeIntegration(msg,Lp,Fp_new,Fe_new,Tstar_v,P,state_new,dt_aim,e,i,g,Temperature,Fg_aim,Fp_current) +! + + if (msg == 'ok') then + subFrac = subFrac + subStep + subStep = min(1.0_pReal-subFrac, subStep*2.0_pReal) ! accelerate + nCutbacks = 0_pInt ! reset cutback counter + success = .true. ! keep current Lp + else + nCutbacks = nCutbacks + 1 ! record additional cutback + maxCutbacks = max(nCutbacks,maxCutbacks)! remember maximum number of cutbacks + subStep = subStep / 2.0_pReal ! cut time step in half + success = .false. ! force Lp interpolation +! if (debugger) then +!$OMP CRITICAL (write2out) + write (6,*) '>>>>>>>>>>>>>>>>>>>> cutback <<<<<<<<<<<<<<<<<<<<<<' + write (6,*) 'Element, Ip:', e, i + write (6,*) msg +!$OMP END CRITICAL (write2out) +! endif +! + endif + enddo ! potential substepping +! +!$OMP CRITICAL (cutback) + debug_cutbackDistribution(min(nCutback,maxCutbacks)+1) = debug_cutbackDistribution(min(nCutback,maxCutbacks)+1)+1 +!$OMP END CRITICAL (cutback) +! +! debugger = .false. + return + end subroutine + +! +!*********************************************************************** +!*** fully-implicit two-level time integration *** +!*** based on a residuum in Lp and intermediate *** +!*** acceleration of the Newton-Raphson correction *** +!*********************************************************************** + SUBROUTINE TimeIntegration(& + msg,& ! return message + Lpguess,& ! guess of plastic velocity gradient + Fp_new,& ! new plastic deformation gradient + Fe_new,& ! new "elastic" deformation gradient + Tstar_v,& ! Stress vector + P,& ! 1nd PK stress (taken as initial guess if /= 0) + state,& ! current microstructure at end of time inc (taken as guess if /= 0) + dt,& ! time increment + cp_en,& ! element number + ip,& ! integration point number + grain,& ! grain number + Temperature,& ! temperature + Fg_new,& ! new total def gradient + Fp_old) ! former plastic def gradient +! state_current) ! former microstructure + use prec + use debug + use mesh, only: mesh_element + use constitutive, only: constitutive_Nstatevars,constitutive_Microstructure,& + constitutive_homogenizedC,constitutive_LpAndItsTangent + use math + use IO + implicit none +! + character(len=*) msg + logical failed + integer(pInt) cp_en, ip, grain + integer(pInt) iInner,dummy, i,j,k,l,m,n + real(pReal) dt, Temperature, det, p_hydro, leapfrog,maxleap + real(pReal), dimension(6) :: Tstar_v + real(pReal), dimension(9,9) :: dLp,dTdLp,dRdLp,invdRdLp,eye2 + real(pReal), dimension(6,6) :: C_66 + real(pReal), dimension(3,3) :: Fg_new,Fp_new,invFp_new,Fp_old,invFp_old,Fe_new + real(pReal), dimension(3,3) :: P !,Tstar + real(pReal), dimension(3,3) :: Lp,Lpguess,Lpguess_old,Rinner,Rinner_old,A,B,BT,AB,BTA + real(pReal), dimension(3,3,3,3) :: C + real(pReal), dimension(constitutive_Nstatevars(grain, ip, cp_en)) :: state +! + msg = 'ok' ! error-free so far + eye2 = math_identity2nd(9) + + call math_invert3x3(Fp_old,invFp_old,det,failed) ! inversion of Fp_old + if (failed) then + msg = 'inversion Fp_old' + return + endif + + A = math_mul33x33(transpose(invFp_old), math_mul33x33(transpose(Fg_new),math_mul33x33(Fg_new,invFp_old))) +! +! if (all(state == 0.0_pReal)) state = state_current ! former state guessed, if none specified +! iOuter = 0_pInt ! outer counter +! + if (debugger) then +!$OMP CRITICAL (write2out) + write (6,'(a,/,3(3(f12.7,x)/))') 'Fg to be calculated',Fg_new +!$OMP END CRITICAL (write2out) + endif +! +!Outer: do ! outer iteration: State +! iOuter = iOuter+1 +! if (debugger) then +!!$OMP CRITICAL (write2out) +! write (6,'(a,i3)') '---outer ',iOuter +! write (6,'(a,/,3(4(f9.3,x)/))') 'state old / MPa',state_old/1e6_pReal +! write (6,'(a,/,3(4(f9.3,x)/))') 'state / MPa',state/1e6_pReal +! write (6,'(a,/,3(3(f12.7,x)/))') 'Lpguess',Lpguess(1:3,:) +!!$OMP END CRITICAL (write2out) +! endif +! +! if (iOuter > nOuter) then +! msg = 'limit Outer iteration' +!!$OMP CRITICAL (out) +! debug_OuterLoopDistribution(nOuter) = debug_OuterLoopDistribution(nOuter)+1 +!!$OMP END CRITICAL (out) +! return +! endif + call constitutive_Microstructure(state,Temperature,grain,ip,cp_en) + C_66 = constitutive_HomogenizedC(state, grain, ip, cp_en) + C = math_Mandel66to3333(C_66) ! 4th rank elasticity tensor +! + iInner = 0_pInt + leapfrog = 1.0_pReal ! correction as suggested by invdRdLp-step + maxleap = 1024.0_pReal ! preassign maximum acceleration level +! + Lpguess_old = Lpguess ! consider present Lpguess good +! +Inner: do ! inner iteration: Lp + iInner = iInner+1 +! if (debugger) then +!!$OMP CRITICAL (write2out) +! write (6,'(a,i3)') 'inner ',iInner +! if (iInner < 3) then +! write (6,'(a,/,3(3(f12.7,x)/))') 'Lpguess',Lpguess(1:3,:) +! endif +!!$OMP END CRITICAL (write2out) +! endif + if (iInner > nInner) then ! too many loops required + Lpguess = Lpguess_old ! do not trust the last update but resort to former one + msg = 'limit Inner iteration' +!$OMP CRITICAL (in) + debug_InnerLoopDistribution(nInner) = debug_InnerLoopDistribution(nInner)+1 +!$OMP END CRITICAL (in) + return + endif +! + B = math_i3 - dt*Lpguess + BT = transpose(B) + AB = math_mul33x33(A,B) + BTA = math_mul33x33(BT,A) + Tstar_v = 0.5_pReal*math_mul66x6(C_66,math_mandel33to6(math_mul33x33(BT,AB)-math_I3)) +! Tstar = math_Mandel6to33(Tstar_v) + p_hydro=(Tstar_v(1)+Tstar_v(2)+Tstar_v(3))/3.0_pReal + forall(i=1:3) Tstar_v(i) = Tstar_v(i)-p_hydro ! subtract hydrostatic pressure + call constitutive_LpAndItsTangent(Lp,dLp, & + Tstar_v,state,Temperature,grain,ip,cp_en) +! + Rinner = Lpguess - Lp ! update current residuum +! + if (.not.(any(Rinner/=Rinner)) .and. & ! exclude any NaN in residuum + ( (maxval(abs(Rinner)) < abstol_Inner) .or. & ! below abs tol .or. + ( any(abs(dt*Lpguess) > relevantStrain) .and. & ! worth checking? .and. + maxval(abs(Rinner/Lpguess),abs(dt*Lpguess) > relevantStrain) < reltol_Inner & ! below rel tol + ) & + ) & + ) & + exit Inner ! convergence +! +! check for acceleration/deceleration in Newton--Raphson correction +! + if (any(Rinner/=Rinner) .and. & ! NaN occured at regular speed + leapfrog == 1.0) then + Lpguess = Lpguess_old ! restore known good guess + msg = 'NaN present' ! croak for cutback + return + + elseif (leapfrog > 1.0_pReal .and. & ! at fast pace ? + (sum(Rinner*Rinner) > sum(Rinner_old*Rinner_old) .or. & ! worse residuum + sum(Rinner*Rinner_old) < 0.0_pReal) .or. & ! residuum changed sign (overshoot) + any(Rinner/=Rinner) ) then ! NaN + maxleap = 0.5_pReal * leapfrog ! limit next acceleration + leapfrog = 1.0_pReal ! grinding halt + + else ! better residuum + dTdLp = 0.0_pReal ! calc dT/dLp + forall (i=1:3,j=1:3,k=1:3,l=1:3,m=1:3,n=1:3) & + dTdLp(3*(i-1)+j,3*(k-1)+l) = dTdLp(3*(i-1)+j,3*(k-1)+l) + & + C(i,j,l,n)*AB(k,n)+C(i,j,m,l)*BTA(m,k) + dTdLp = -0.5_pReal*dt*dTdLp + dRdLp = eye2 - math_mul99x99(dLp,dTdLp) ! calc dR/dLp + invdRdLp = 0.0_pReal + call math_invert(9,dRdLp,invdRdLp,dummy,failed) ! invert dR/dLp --> dLp/dR + if (failed) then + msg = 'inversion dR/dLp' + if (debugger) then +!$OMP CRITICAL (write2out) + write (6,*) msg + write (6,'(a,/,9(9(e9.3,x)/))') 'dRdLp', dRdLp(1:9,:) + write (6,'(a,/,3(4(f9.3,x)/))') 'state / MPa',state/1e6_pReal + write (6,'(a,/,3(3(f12.7,x)/))') 'Lpguess',Lpguess(1:3,:) + write (6,'(a,/,3(3(e12.7,x)/))') 'Lp',Lp(1:3,:) + write (6,'(a,/,6(f9.3,x))') 'Tstar / MPa',Tstar_v/1e6_pReal +!$OMP END CRITICAL (write2out) + endif + return + endif +! + Rinner_old = Rinner ! remember current residuum + Lpguess_old = Lpguess ! remember current Lp guess + if (iInner > 1 .and. leapfrog < maxleap) & + leapfrog = 2.0_pReal * leapfrog ! accelerate if ok + endif +! + Lpguess = Lpguess_old ! start from current guess + Rinner = Rinner_old ! use current residuum + forall (i=1:3,j=1:3,k=1:3,l=1:3) & ! leapfrog to updated Lpguess + Lpguess(i,j) = Lpguess(i,j) - leapfrog*invdRdLp(3*(i-1)+j,3*(k-1)+l)*Rinner(k,l) + enddo Inner +! +!$OMP CRITICAL (in) + debug_InnerLoopDistribution(iInner) = debug_InnerLoopDistribution(iInner)+1 +!$OMP END CRITICAL (in) +! ROuter = state - state_old - & +! dt*constitutive_dotState(Tstar_v,state,Temperature,& +! grain,ip,cp_en) ! residuum from evolution of microstructure +! state = state - ROuter ! update of microstructure +! +! if (iOuter==nOuter) then +!!$OMP CRITICAL (write2out) +! write (6,*) 'Terminated outer loop at el,ip,grain',cp_en,ip,grain +!!$OMP END CRITICAL (write2out) +! exit Outer +! endif +! if (maxval(abs(Router/state),state /= 0.0_pReal) < reltol_Outer) exit Outer +! enddo Outer +! +!!$OMP CRITICAL (out) +! debug_OuterLoopDistribution(iOuter) = debug_OuterLoopDistribution(iOuter)+1 +!!$OMP END CRITICAL (out) + invFp_new = math_mul33x33(invFp_old,B) + call math_invert3x3(invFp_new,Fp_new,det,failed) + if (failed) then + msg = 'inversion Fp_new^-1' + return + endif +! +! if (wantsConstitutiveResults) then ! get the post_results upon request +! results = 0.0_pReal +! results = constitutive_post_results(Tstar_v,state,Temperature,dt,grain,ip,cp_en) +! endif +! + Fp_new = Fp_new*det**(1.0_pReal/3.0_pReal) ! regularize Fp by det = det(InvFp_new) !! + forall (i=1:3) Tstar_v(i) = Tstar_v(i)+p_hydro ! add hydrostatic component back + Fe_new = math_mul33x33(Fg_new,invFp_new) ! calc resulting Fe +! P = math_mul33x33(Fe_new,math_mul33x33(Tstar,transpose(invFp_new))) ! first PK stress + P = math_mul33x33(Fe_new,math_mul33x33(math_Mandel6to33(Tstar_v),transpose(invFp_new))) ! first PK stress + + return +! + END SUBROUTINE +! + END MODULE +!############################################################## + diff --git a/trunk/crystallite.f90 b/trunk/crystallite.f90 deleted file mode 100644 index 5e80bfb29..000000000 --- a/trunk/crystallite.f90 +++ /dev/null @@ -1,398 +0,0 @@ -!############################################################## - MODULE crystallite -!############################################################## -! *** Solution at single crystallite level *** -! -CONTAINS -! -! -!******************************************************************** -! Calculates the stress for a single component -!******************************************************************** -!*********************************************************************** -!*** calculation of stress (P), stiffness (dPdF), *** -!*** and announcment of any *** -!*** acceleration of the Newton-Raphson correction *** -!*********************************************************************** - subroutine SingleCrystallite(& - msg,& ! return message - P,& ! first PK stress - dPdF,& ! consistent tangent - post_results,& ! plot results from constitutive model - Lp,& ! plastic velocity gradient - Fp_new,& ! new plastic deformation gradient - Fe_new,& ! new "elastic" deformation gradient - state_new,& ! new state variable array - dt,& ! time increment - cp_en,& ! element number - ip,& ! integration point number - grain,& ! grain number - updateJaco,& ! update of Jacobian required - Temperature,& ! temperature of crystallite - Fg_new,& ! new global deformation gradient - Fg_old,& ! old global deformation gradient - Fp_old,& ! old plastic deformation gradient - state_old) ! old state variable array -! - use prec, only: pReal,pInt,pert_Fg,subStepMin, nCutback - use debug - use constitutive, only: constitutive_Nstatevars,constitutive_Nresults - use mesh, only: mesh_element - use math - - implicit none -! - character(len=*) msg - logical updateJaco,error,success,guessNew - integer(pInt) cp_en,ip,grain,k,l, nCutbacks, maxCutbacks - real(pReal) Temperature - real(pReal) dt,dt_aim,subFrac,subStep,det - real(pReal), dimension(3,3) :: Lp,Lp_interpolated,Lp_pert,inv - real(pReal), dimension(3,3) :: Fg_old,Fg_current,Fg_new,Fg_pert,Fg_aim,deltaFg - real(pReal), dimension(3,3) :: Fp_old,Fp_current,Fp_new,Fp_pert - real(pReal), dimension(3,3) :: Fe_current,Fe_new,Fe_pert - real(pReal), dimension(3,3) :: P,P_pert - real(pReal), dimension(3,3,3,3) :: dPdF - real(pReal), dimension(constitutive_Nstatevars(grain,ip,cp_en)) :: state_old,state_new - real(pReal), dimension(constitutive_Nstatevars(grain,ip,cp_en)) :: state_current,state_bestguess,state_pert - real(pReal), dimension(constitutive_Nresults(grain,ip,cp_en)) :: post_results - - deltaFg = Fg_new - Fg_old - subFrac = 0.0_pReal - subStep = 1.0_pReal - nCutbacks = 0_pInt - maxCutbacks = 0_pInt - Fg_current = Fg_old ! initialize to start of inc - Fp_current = Fp_old - call math_invert3x3(Fp_old,inv,det,error) - Fe_current = math_mul33x33(Fg_old,inv) - state_current = state_old - success = .false. ! pretend cutback - dt_aim = 0.0_pReal ! prevent initial Lp interpolation - -! begin the cutback loop - do while (subStep > subStepMin) ! continue until finished or too much cut backing - if (success) then ! wind forward - Fg_current = Fg_aim - Fe_current = Fe_new - Fp_current = Fp_new - state_current = state_new - elseif (dt_aim > 0.0_pReal) then - call math_invert3x3(Fg_aim,inv,det,error) ! inv of Fg_aim - Lp_interpolated = 0.5_pReal*Lp + & - 0.5_pReal*(math_I3 - math_mul33x33(Fp_current,& - math_mul33x33(inv,Fe_current)))/dt_aim ! interpolate Lp and L - if (debugger) then -!$OMP CRITICAL (write2out) - write (6,*) 'Lp interpolation' - write (6,'(a,/,3(3(f12.7,x)/))') 'from',Lp(1:3,:) - write (6,'(a,/,3(3(f12.7,x)/))') 'to',Lp_interpolated(1:3,:) -!$OMP END CRITICAL (write2out) - endif - Lp = Lp_interpolated - endif - - Fg_aim = Fg_current + subStep*deltaFg ! aim for Fg - dt_aim = subStep*dt ! aim for dt - if (debugger) then -!$OMP CRITICAL (write2out) - write (6,*) 'using these values' - write (6,'(a,/,3(4(f9.3,x)/))') 'state current / MPa',state_current/1e6_pReal - write (6,'(a,/,3(4(f9.3,x)/))') 'state new / MPa',state_new/1e6_pReal - write (6,'(a,/,3(3(f12.7,x)/))') 'Fe current',Fe_current(1:3,:) - write (6,'(a,/,3(3(f12.7,x)/))') 'Fp current',Fp_current(1:3,:) - write (6,'(a,/,3(3(f12.7,x)/))') 'Lp (old=new guess)',Lp(1:3,:) - write (6,'(a20,f,x,a2,x,f)') 'integrating from ',subFrac,'to',(subFrac+subStep) -!$OMP END CRITICAL (write2out) - endif - call TimeIntegration(msg,Lp,Fp_new,Fe_new,P,state_new,post_results, & ! def gradients and PK2 at end of time step - maxval(abs(Fg_aim-Fg_new)) < relevantStrain, & ! post results only if asking for final values - dt_aim,cp_en,ip,grain,Temperature,Fg_aim,Fp_current,state_current) - - if (msg == 'ok') then - subFrac = subFrac + subStep - subStep = min(1.0_pReal-subFrac, subStep*2.0_pReal) ! accelerate - nCutbacks = 0_pInt ! reset cutback counter - success = .true. ! keep current Lp - else - nCutbacks = nCutbacks + 1 ! record additional cutback - maxCutbacks = max(nCutbacks,maxCutbacks)! remember maximum number of cutbacks - subStep = subStep / 2.0_pReal ! cut time step in half - success = .false. ! force Lp interpolation - if (debugger) then -!$OMP CRITICAL (write2out) - write (6,*) '>>>>>>>>>>>>>>>>>>>> cutback <<<<<<<<<<<<<<<<<<<<<<' -!$OMP END CRITICAL (write2out) - endif - - endif - enddo ! potential substepping -! -!$OMP CRITICAL (cutback) - debug_cutbackDistribution(min(nCutback,maxCutbacks)+1) = debug_cutbackDistribution(min(nCutback,maxCutbacks)+1)+1 -!$OMP END CRITICAL (cutback) -! - if (msg /= 'ok') return ! solution not reached --> report back - if (updateJaco) then ! consistent tangent using - if (debugger) then -!$OMP CRITICAL (write2out) - write (6,*) 'Jacobian calc' -!$OMP END CRITICAL (write2out) - endif - do k=1,3 - do l=1,3 - Fg_pert = Fg_new ! initialize perturbed Fg - Fg_pert(k,l) = Fg_pert(k,l) + pert_Fg ! perturb single component - Lp_pert = Lp - state_pert = state_new ! initial guess from end of time step - call TimeIntegration(msg,Lp_pert,Fp_pert,Fe_pert,P_pert,state_pert,post_results,.false., & ! def gradients and PK2 at end of time step - dt_aim,cp_en,ip,grain,Temperature,Fg_pert,Fp_current,state_current) - - if (msg == 'ok') & - dPdF(:,:,k,l) = (P_pert-P)/pert_Fg ! constructing tangent dP_ij/dFg_kl only if valid forward difference - ! otherwise leave component unchanged - enddo - enddo - endif -! - - msg = 'ok' ! a new consistent tangent was computed even if msg was not ok for all components - -! - return -! - END SUBROUTINE -! -!*********************************************************************** -!*** fully-implicit two-level time integration *** -!*** based on a residuum in Lp and intermediate *** -!*** acceleration of the Newton-Raphson correction *** -!*********************************************************************** - SUBROUTINE TimeIntegration(& - msg,& ! return message - Lpguess,& ! guess of plastic velocity gradient - Fp_new,& ! new plastic deformation gradient - Fe_new,& ! new "elastic" deformation gradient - P,& ! 1nd PK stress (taken as initial guess if /= 0) - state,& ! current microstructure at end of time inc (taken as guess if /= 0) - results,& ! post results from constitutive - wantsConstitutiveResults,& ! its flag -! - dt,& ! time increment - cp_en,& ! element number - ip,& ! integration point number - grain,& ! grain number - Temperature,& ! temperature - Fg_new,& ! new total def gradient - Fp_old,& ! former plastic def gradient - state_old) ! former microstructure - use prec - use debug - use mesh, only: mesh_element - use constitutive, only: constitutive_Nstatevars,& - constitutive_homogenizedC,constitutive_dotState,constitutive_LpAndItsTangent,& - constitutive_Nresults,constitutive_Microstructure,constitutive_post_results - use math - use IO - implicit none -! - character(len=*) msg - logical failed,wantsConstitutiveResults - integer(pInt) cp_en, ip, grain - integer(pInt) iOuter,iInner,dummy, i,j,k,l,m,n - real(pReal) dt, Temperature, det, p_hydro, leapfrog,maxleap - real(pReal), dimension(6) :: Tstar_v - real(pReal), dimension(9,9) :: dLp,dTdLp,dRdLp,invdRdLp,eye2 - real(pReal), dimension(6,6) :: C_66 - real(pReal), dimension(3,3) :: Fg_new,Fp_new,invFp_new,Fp_old,invFp_old,Fe_new - real(pReal), dimension(3,3) :: P,Tstar - real(pReal), dimension(3,3) :: Lp,Lpguess,Lpguess_old,Rinner,Rinner_old,A,B,BT,AB,BTA - real(pReal), dimension(3,3,3,3) :: C - real(pReal), dimension(constitutive_Nstatevars(grain, ip, cp_en)) :: state_old,state,ROuter - real(pReal), dimension(constitutive_Nresults(grain,ip,cp_en)) :: results -! - msg = 'ok' ! error-free so far - eye2 = math_identity2nd(9) - - call math_invert3x3(Fp_old,invFp_old,det,failed) ! inversion of Fp_old - if (failed) then - msg = 'inversion Fp_old' - return - endif - - A = math_mul33x33(transpose(invFp_old), math_mul33x33(transpose(Fg_new),math_mul33x33(Fg_new,invFp_old))) -! - if (all(state == 0.0_pReal)) state = state_old ! former state guessed, if none specified - iOuter = 0_pInt ! outer counter -! - if (debugger) then -!$OMP CRITICAL (write2out) - write (6,'(a,/,3(3(f12.7,x)/))') 'Fg to be calculated',Fg_new -!$OMP END CRITICAL (write2out) - endif -! -Outer: do ! outer iteration: State - iOuter = iOuter+1 - if (debugger) then -!$OMP CRITICAL (write2out) - write (6,'(a,i3)') '---outer ',iOuter - write (6,'(a,/,3(4(f9.3,x)/))') 'state old / MPa',state_old/1e6_pReal - write (6,'(a,/,3(4(f9.3,x)/))') 'state / MPa',state/1e6_pReal - write (6,'(a,/,3(3(f12.7,x)/))') 'Lpguess',Lpguess(1:3,:) -!$OMP END CRITICAL (write2out) - endif - - if (iOuter > nOuter) then - msg = 'limit Outer iteration' -!$OMP CRITICAL (out) - debug_OuterLoopDistribution(nOuter) = debug_OuterLoopDistribution(nOuter)+1 -!$OMP END CRITICAL (out) - return - endif - call constitutive_Microstructure(state,Temperature,grain,ip,cp_en) - C_66 = constitutive_HomogenizedC(state, grain, ip, cp_en) - C = math_Mandel66to3333(C_66) ! 4th rank elasticity tensor -! - iInner = 0_pInt - leapfrog = 1.0_pReal ! correction as suggested by invdRdLp-step - maxleap = 1024.0_pReal ! preassign maximum acceleration level -! - Lpguess_old = Lpguess ! consider present Lpguess good -! -Inner: do ! inner iteration: Lp - iInner = iInner+1 - if (debugger) then -!$OMP CRITICAL (write2out) - write (6,'(a,i3)') 'inner ',iInner - if (wantsConstitutiveResults .and. iOuter == 1 .and. iInner < 3) then - write (6,'(a,/,3(3(f12.7,x)/))') 'Lpguess',Lpguess(1:3,:) - endif -!$OMP END CRITICAL (write2out) - endif - if (iInner > nInner) then ! too many loops required - Lpguess = Lpguess_old ! do not trust the last update but resort to former one - msg = 'limit Inner iteration' -!$OMP CRITICAL (in) - debug_InnerLoopDistribution(nInner) = debug_InnerLoopDistribution(nInner)+1 -!$OMP END CRITICAL (in) - return - endif -! - B = math_i3 - dt*Lpguess - BT = transpose(B) - AB = math_mul33x33(A,B) - BTA = math_mul33x33(BT,A) - Tstar_v = 0.5_pReal*math_mul66x6(C_66,math_mandel33to6(math_mul33x33(BT,AB)-math_I3)) - Tstar = math_Mandel6to33(Tstar_v) - p_hydro=(Tstar_v(1)+Tstar_v(2)+Tstar_v(3))/3.0_pReal - forall(i=1:3) Tstar_v(i) = Tstar_v(i)-p_hydro ! subtract hydrostatic pressure - call constitutive_LpAndItsTangent(Lp,dLp, & - Tstar_v,state,Temperature,grain,ip,cp_en) -! - Rinner = Lpguess - Lp ! update current residuum -! - if (.not.(any(Rinner/=Rinner)) .and. & ! exclude any NaN in residuum - ( (maxval(abs(Rinner)) < abstol_Inner) .or. & ! below abs tol .or. - ( any(abs(dt*Lpguess) > relevantStrain) .and. & ! worth checking? .and. - maxval(abs(Rinner/Lpguess),abs(dt*Lpguess) > relevantStrain) < reltol_Inner & ! below rel tol - ) & - ) & - ) & - exit Inner ! convergence -! -! check for acceleration/deceleration in Newton--Raphson correction -! - if (any(Rinner/=Rinner) .and. & ! NaN occured at regular speed - leapfrog == 1.0) then - Lpguess = Lpguess_old ! restore known good guess - msg = 'NaN present' ! croak for cutback - return - - elseif (leapfrog > 1.0_pReal .and. & ! at fast pace ? - (sum(Rinner*Rinner) > sum(Rinner_old*Rinner_old) .or. & ! worse residuum - sum(Rinner*Rinner_old) < 0.0_pReal) .or. & ! residuum changed sign (overshoot) - any(Rinner/=Rinner) ) then ! NaN - maxleap = 0.5_pReal * leapfrog ! limit next acceleration - leapfrog = 1.0_pReal ! grinding halt - - else ! better residuum - dTdLp = 0.0_pReal ! calc dT/dLp - forall (i=1:3,j=1:3,k=1:3,l=1:3,m=1:3,n=1:3) & - dTdLp(3*(i-1)+j,3*(k-1)+l) = dTdLp(3*(i-1)+j,3*(k-1)+l) + & - C(i,j,l,n)*AB(k,n)+C(i,j,m,l)*BTA(m,k) - dTdLp = -0.5_pReal*dt*dTdLp - dRdLp = eye2 - math_mul99x99(dLp,dTdLp) ! calc dR/dLp - invdRdLp = 0.0_pReal - call math_invert(9,dRdLp,invdRdLp,dummy,failed) ! invert dR/dLp --> dLp/dR - if (failed) then - msg = 'inversion dR/dLp' - if (debugger) then -!$OMP CRITICAL (write2out) - write (6,*) msg - write (6,'(a,/,9(9(e9.3,x)/))') 'dRdLp', dRdLp(1:9,:) - write (6,'(a,/,3(4(f9.3,x)/))') 'state / MPa',state/1e6_pReal - write (6,'(a,/,3(3(f12.7,x)/))') 'Lpguess',Lpguess(1:3,:) - write (6,'(a,/,3(3(e12.7,x)/))') 'Lp',Lp(1:3,:) - write (6,'(a,/,6(f9.3,x))') 'Tstar / MPa',Tstar_v/1e6_pReal -!$OMP END CRITICAL (write2out) - endif - return - endif -! - Rinner_old = Rinner ! remember current residuum - Lpguess_old = Lpguess ! remember current Lp guess - if (iInner > 1 .and. leapfrog < maxleap) & - leapfrog = 2.0_pReal * leapfrog ! accelerate if ok - endif -! - Lpguess = Lpguess_old ! start from current guess - Rinner = Rinner_old ! use current residuum - forall (i=1:3,j=1:3,k=1:3,l=1:3) & ! leapfrog to updated Lpguess - Lpguess(i,j) = Lpguess(i,j) - leapfrog*invdRdLp(3*(i-1)+j,3*(k-1)+l)*Rinner(k,l) - enddo Inner -! -!$OMP CRITICAL (in) - debug_InnerLoopDistribution(iInner) = debug_InnerLoopDistribution(iInner)+1 -!$OMP END CRITICAL (in) - ROuter = state - state_old - & - dt*constitutive_dotState(Tstar_v,state,Temperature,& - grain,ip,cp_en) ! residuum from evolution of microstructure - state = state - ROuter ! update of microstructure - - if (iOuter==nOuter) then -!!$OMP CRITICAL (write2out) - write (6,*) 'Terminated outer loop at el,ip,grain',cp_en,ip,grain -!!$OMP END CRITICAL (write2out) - exit Outer - endif - if (maxval(abs(Router/state),state /= 0.0_pReal) < reltol_Outer) exit Outer - enddo Outer -! -!$OMP CRITICAL (out) - debug_OuterLoopDistribution(iOuter) = debug_OuterLoopDistribution(iOuter)+1 -!$OMP END CRITICAL (out) - invFp_new = math_mul33x33(invFp_old,B) - call math_invert3x3(invFp_new,Fp_new,det,failed) - if (failed) then - msg = 'inversion Fp_new^-1' - return - endif -! - if (wantsConstitutiveResults) then ! get the post_results upon request - results = 0.0_pReal - results = constitutive_post_results(Tstar_v,state,Temperature,dt,grain,ip,cp_en) - endif -! - Fp_new = Fp_new*det**(1.0_pReal/3.0_pReal) ! regularize Fp by det = det(InvFp_new) !! - forall (i=1:3) Tstar_v(i) = Tstar_v(i)+p_hydro ! add hydrostatic component back - Fe_new = math_mul33x33(Fg_new,invFp_new) ! calc resulting Fe - P = math_mul33x33(Fe_new,math_mul33x33(Tstar,transpose(invFp_new))) ! first PK stress - - return -! - END SUBROUTINE -! -! - END MODULE -!############################################################## - diff --git a/trunk/mpie_cpfem_marc2005r3.f90 b/trunk/mpie_cpfem_marc2005r3.f90 deleted file mode 100644 index 6514cd0b7..000000000 --- a/trunk/mpie_cpfem_marc2005r3.f90 +++ /dev/null @@ -1,332 +0,0 @@ -!******************************************************************** -! Material subroutine for MSC.Marc Version 0.1 -! -! written by F. Roters, P. Eisenlohr, L. Hantcherli, W.A. Counts -! MPI fuer Eisenforschung, Duesseldorf -! -! last modified: 08.11.2007 -!******************************************************************** -! Usage: -! - choose material as hypela2 -! - set statevariable 2 to index of material -! - set statevariable 3 to index of texture -! - choose output of user variables if desired -! - make sure the file "mattex.mpie" exists in the working -! directory -! - use nonsymmetric option for solver (e.g. direct -! profile or multifrontal sparse, the latter seems -! to be faster!) -!******************************************************************** -! Marc subroutines used: -! - hypela2 -! - plotv -! - quit -!******************************************************************** -! Marc common blocks included: -! - concom: lovl, ncycle, inc, incsub -! - creeps: timinc -!******************************************************************** -! - include "prec.f90" ! uses nothing else - include "debug.f90" ! uses prec - include "math.f90" ! uses prec - include "IO.f90" ! uses prec, debug, math - include "FEsolving.f90" ! uses prec, IO - include "mesh.f90" ! uses prec, IO, math, FEsolving - include "lattice.f90" ! uses prec, math - include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug - include "crystallite.f90" ! uses prec, debug, constitutive, mesh, math, IO - include "CPFEM.f90" ! uses prec, math, mesh, constitutive, FEsolving, debug, lattice, IO, crystallite -! - - SUBROUTINE hypela2(d,g,e,de,s,t,dt,ngens,n,nn,kcus,matus,ndi,& - nshear,disp,dispt,coord,ffn,frotn,strechn,eigvn,ffn1,& - frotn1,strechn1,eigvn1,ncrd,itel,ndeg,ndm,& - nnode,jtype,lclass,ifr,ifu) -!******************************************************************** -! This is the Marc material routine -!******************************************************************** -! -! ************* user subroutine for defining material behavior ************** -! -! -! CAUTION : Due to calculation of the Deformation gradients, Stretch Tensors and -! Rotation tensors at previous and current states, the analysis can be -! computationally expensive. Please use the user subroutine -> hypela -! if these kinematic quantities are not needed in the constitutive model -! -! -! IMPORTANT NOTES : -! -! (1) F,R,U are only available for continuum and membrane elements (not for -! shells and beams). -! -! (2) For total Lagrangian formulation use the -> 'Elasticity,1' card(= -! total Lagrange with large disp) in the parameter section of input deck. -! For updated Lagrangian formulation use the -> 'Plasticity,3' card(= -! update+finite+large disp+constant d) in the parameter section of -! input deck. -! -! -! d stress strain law to be formed -! g change in stress due to temperature effects -! e total elastic strain -! de increment of strain -! s stress - should be updated by user -! t state variables (comes in at t=n, must be updated -! to have state variables at t=n+1) -! dt ! Marc common blocks are in fixed format so they have to be pasted in here -! -! Marc common blocks are in fixed format so they have to be pasted in here -! Beware of changes in newer Marc versions -- these are from 2005r3 -! concom is needed for inc, subinc, ncycle, lovl -! include 'concom' - common/concom/ & - iacous, iasmbl, iautth, ibear, icompl, iconj, icreep, ideva(50), idyn, idynt,& - ielas, ielcma, ielect, iform, ifour, iharm, ihcps, iheat, iheatt, ihresp,& - ijoule, ilem, ilnmom, iloren, inc, incext, incsub, ipass, iplres, ipois,& - ipoist, irpflo, ismall, ismalt, isoil, ispect, ispnow, istore, iswep, ithcrp,& - itherm, iupblg, iupdat, jacflg, jel, jparks, largst, lfond, loadup, loaduq,& - lodcor, lovl, lsub, magnet, ncycle, newtnt, newton, noshr, linear, ivscpl,& - icrpim, iradrt, ipshft, itshr, iangin, iupmdr, iconjf, jincfl, jpermg, jhour,& - isolvr, jritz, jtable, jshell, jdoubl, jform, jcentr, imini, kautth, iautof,& - ibukty, iassum, icnstd, icnstt, kmakmas, imethvp,iradrte,iradrtp, iupdate,iupdatp,& - ncycnt, marmen ,idynme, ihavca, ispf, kmini, imixed, largtt, kdoela, iautofg,& - ipshftp,idntrc, ipore, jtablm, jtablc, isnecma,itrnspo,imsdif, jtrnspo,mcnear,& - imech, imecht, ielcmat, ielectt,magnett, imsdift,noplas, jtabls, jactch, jtablth,& - kgmsto ,jpzo, ifricsh, iremkin,iremfor, ishearp,jspf, machining, jlshell,icompsol,& - iupblgfo,jcondir,nstcrp, nactive,ipassref, nstspnt,ibeart,icheckmpc, noline, icuring,& - ishrink,ioffsflg,isetoff, iharmt, inc_incdat, iautspc,ibrake - -! creeps is needed for timinc (time increment) -! include 'creeps' - common/marc_creeps/ & - cptim,timinc,timinc_p,timinc_s,timincm,timinc_a,timinc_b,creept(33),icptim,icfte,icfst,& - icfeq,icftm,icetem,mcreep,jcreep,icpa,icftmp,icfstr,icfqcp,icfcpm,icrppr,icrcha,icpb,iicpmt,iicpa - - - increment of state variables -! ngens size of stress - strain law -! n element number -! nn integration point number -! kcus(1) layer number -! kcus(2) internal layer number -! matus(1) user material identification number -! matus(2) internal material identification number -! ndi number of direct components -! nshear number of shear components -! disp incremental displacements -! dispt displacements at t=n (at assembly, lovl=4) and -! displacements at t=n+1 (at stress recovery, lovl=6) -! coord coordinates -! ncrd number of coordinates -! ndeg number of degrees of freedom -! itel dimension of F and R, either 2 or 3 -! nnode number of nodes per element -! jtype element type -! lclass element class -! ifr set to 1 if R has been calculated -! ifu set to 1 if strech has been calculated -! -! at t=n : -! -! ffn deformation gradient -! frotn rotation tensor -! strechn square of principal stretch ratios, lambda(i) -! eigvn(i,j) i principal direction components for j eigenvalues -! -! at t=n+1 : -! -! ffn1 deformation gradient -! frotn1 rotation tensor -! strechn1 square of principal stretch ratios, lambda(i) -! eigvn1(i,j) i principal direction components for j eigenvalues -! -! The following operation obtains U (stretch tensor) at t=n+1 : -! -! call scla(un1,0.d0,itel,itel,1) -! do 3 k=1,3 -! do 2 i=1,3 -! do 1 j=1,3 -! un1(i,j)=un1(i,j)+dsqrt(strechn1(k))*eigvn1(i,k)*eigvn1(j,k) -!1 continue -!2 continue -!3 continue -! - - use prec, only: pReal,pInt, ijaco - use FEsolving - use CPFEM, only: CPFEM_general - use math, only: invnrmMandel - - implicit real(pReal) (a-h,o-z) - - integer(pInt) computationMode - - - - dimension e(*),de(*),t(*),dt(*),g(*),d(ngens,*),s(*), n(2),coord(ncrd,*),disp(ndeg,*),matus(2),dispt(ndeg,*),ffn(itel,*),& - frotn(itel,*),strechn(itel),eigvn(itel,*),ffn1(itel,*),frotn1(itel,*),strechn1(itel),eigvn1(itel,*),kcus(2) - - - -! Marc common blocks are in fixed format so they have to be pasted in here -! Beware of changes in newer Marc versions -- these are from 2005r3 -! concom is needed for inc, subinc, ncycle, lovl -! include 'concom' - common/marc_concom/ & - iacous, iasmbl, iautth, ibear, icompl, iconj, icreep, ideva(50), idyn, idynt,& - ielas, ielcma, ielect, iform, ifour, iharm, ihcps, iheat, iheatt, ihresp,& - ijoule, ilem, ilnmom, iloren, inc, incext, incsub, ipass, iplres, ipois,& - ipoist, irpflo, ismall, ismalt, isoil, ispect, ispnow, istore, iswep, ithcrp,& - itherm, iupblg, iupdat, jacflg, jel, jparks, largst, lfond, loadup, loaduq,& - lodcor, lovl, lsub, magnet, ncycle, newtnt, newton, noshr, linear, ivscpl,& - icrpim, iradrt, ipshft, itshr, iangin, iupmdr, iconjf, jincfl, jpermg, jhour,& - isolvr, jritz, jtable, jshell, jdoubl, jform, jcentr, imini, kautth, iautof,& - ibukty, iassum, icnstd, icnstt, kmakmas, imethvp,iradrte,iradrtp, iupdate,iupdatp,& - ncycnt, marmen ,idynme, ihavca, ispf, kmini, imixed, largtt, kdoela, iautofg,& - ipshftp,idntrc, ipore, jtablm, jtablc, isnecma,itrnspo,imsdif, jtrnspo,mcnear,& - imech, imecht, ielcmat, ielectt,magnett, imsdift,noplas, jtabls, jactch, jtablth,& - kgmsto ,jpzo, ifricsh, iremkin,iremfor, ishearp,jspf, machining, jlshell,icompsol,& - iupblgfo,jcondir,nstcrp, nactive,ipassref, nstspnt,ibeart,icheckmpc, noline, icuring,& - ishrink,ioffsflg,isetoff, ioffsetm,iharmt, inc_incdat,iautspc,ibrake, icbush ,istream_input,& - iprsinp,ivlsinp,ifirst_time,ipin_m,jgnstr_glb, imarc_return,iqvcinp,nqvceid,istpnx,imicro1 - -! creeps is needed for timinc (time increment) -! include 'creeps' - common/marc_creeps/ & - cptim,timinc,timinc_p,timinc_s,timincm,timinc_a,timinc_b,creept(33),icptim,icfte,icfst,& - icfeq,icftm,icetem,mcreep,jcreep,icpa,icftmp,icfstr,icfqcp,icfcpm,icrppr,icrcha,icpb,iicpmt,iicpa - - - if (inc == 0) then - cycleCounter = 0 - else - if (theCycle > ncycle) cycleCounter = 0 ! reset counter for each cutback - if (theCycle /= ncycle .or. theLovl /= lovl) cycleCounter = cycleCounter+1 ! ping pong - endif - if (cptim > theTime .or. theInc /= inc) then ! reached convergence - lastIncConverged = .true. - outdatedByNewInc = .true. - endif - - if (mod(cycleCounter,2) /= 0) computationMode = 4 ! recycle - if (mod(cycleCounter,4) == 2) computationMode = 3 ! collect - if (mod(cycleCounter,4) == 0) computationMode = 2 ! compute - if (computationMode == 4 .and. ncycle == 0 .and. .not. lastIncConverged) & - computationMode = 6 ! recycle but restore known good consistent tangent - if (computationMode == 4 .and. lastIncConverged) then - computationMode = 5 ! recycle and record former consistent tangent - lastIncConverged = .false. - endif - if (computationMode == 2 .and. outdatedByNewInc) then - computationMode = 1 ! compute and age former results - outdatedByNewInc = .false. - endif - - theTime = cptim ! record current starting time - theInc = inc ! record current increment number - theCycle = ncycle ! record current cycle count - theLovl = lovl ! record current lovl - - call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,n(1),nn,s,mod(theCycle,2_pInt*ijaco)==0,d,ngens) - - -! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13 -! Marc: 11, 22, 33, 12, 23, 13 - forall(i=1:ngens) d(1:ngens,i) = invnrmMandel(i)*d(1:ngens,i)*invnrmMandel(1:ngens) - s(1:ngens) = s(1:ngens)*invnrmMandel(1:ngens) - if(symmetricSolver) d(1:ngens,1:ngens) = 0.5_pReal*(d(1:ngens,1:ngens)+transpose(d(1:ngens,1:ngens))) - return - - END SUBROUTINE -! - - SUBROUTINE plotv(v,s,sp,etot,eplas,ecreep,t,m,nn,layer,ndi,nshear,jpltcd) -!******************************************************************** -! This routine sets user defined output variables for Marc -!******************************************************************** -! -! select a variable contour plotting (user subroutine). -! -! v variable -! s (idss) stress array -! sp stresses in preferred direction -! etot total strain (generalized) -! eplas total plastic strain -! ecreep total creep strain -! t current temperature -! m element number -! nn integration point number -! layer layer number -! ndi (3) number of direct stress components -! nshear (3) number of shear stress components -! -!******************************************************************** - use prec, only: pReal,pInt - use CPFEM, only: CPFEM_results, CPFEM_Nresults - use constitutive, only: constitutive_maxNresults - use mesh, only: mesh_FEasCP - implicit none -! - real(pReal) s(*),etot(*),eplas(*),ecreep(*),sp(*) - real(pReal) v, t(*) - integer(pInt) m, nn, layer, ndi, nshear, jpltcd -! -! assign result variable - v=CPFEM_results(mod(jpltcd-1_pInt, CPFEM_Nresults+constitutive_maxNresults)+1_pInt,& - (jpltcd-1_pInt)/(CPFEM_Nresults+constitutive_maxNresults)+1_pInt,& - nn, mesh_FEasCP('elem', m)) - return - END SUBROUTINE -! -! -! subroutine utimestep(timestep,timestepold,icall,time,timeloadcase) -!******************************************************************** -! This routine modifies the addaptive time step of Marc -!******************************************************************** -! use prec, only: pReal,pInt -! use CPFEM, only : CPFEM_timefactor_max -! implicit none -! -! real(pReal) timestep, timestepold, time,timeloadcase -! integer(pInt) icall -! -! user subroutine for modifying the time step in auto step -! -! timestep : the current time step as suggested by marc -! to be modified in this routine -! timestepold : the current time step before it was modified by marc -! icall : =1 for setting the initial time step -! =2 if this routine is called during an increment -! =3 if this routine is called at the beginning -! of the increment -! time : time at the start of the current increment -! timeloadcase: time period of the current load case -! -! it is in general not recommended to increase the time step -! during the increment. -! this routine is called right after the time step has (possibly) -! been updated by marc. -! -! user coding -! reduce timestep during increment in case mpie_timefactor is too large -! if(icall==2_pInt) then -! if(mpie_timefactor_max>1.25_pReal) then -! timestep=min(timestep,timestepold*0.8_pReal) -! end if -! return -! modify timestep at beginning of new increment -! else if(icall==3_pInt) then -! if(mpie_timefactor_max<=0.8_pReal) then -! timestep=min(timestep,timestepold*1.25_pReal) -! else if (mpie_timefactor_max<=1.0_pReal) then -! timestep=min(timestep,timestepold/mpie_timefactor_max) -! else if (mpie_timefactor_max<=1.25_pReal) then -! timestep=min(timestep,timestepold*1.01_pReal) -! else -! timestep=min(timestep,timestepold*0.8_pReal) -! end if -! end if -! return -! end diff --git a/trunk/mpie_cpfem_marc2007r1.f90 b/trunk/mpie_cpfem_marc2007r1.f90 index 89e4871d0..465be638f 100644 --- a/trunk/mpie_cpfem_marc2007r1.f90 +++ b/trunk/mpie_cpfem_marc2007r1.f90 @@ -4,7 +4,7 @@ ! written by F. Roters, P. Eisenlohr, L. Hantcherli, W.A. Counts ! MPI fuer Eisenforschung, Duesseldorf ! -! last modified: 09.07.2008 +! last modified: 27.11.2008 !******************************************************************** ! Usage: ! - choose material as hypela2 @@ -35,7 +35,7 @@ include "mesh.f90" ! uses prec, IO, math, FEsolving include "lattice.f90" ! uses prec, math include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug - include "crystallite.f90" ! uses prec, debug, constitutive, mesh, math, IO +! include "crystallite.f90" ! uses prec, debug, constitutive, mesh, math, IO include "CPFEM.f90" ! uses prec, math, mesh, constitutive, FEsolving, debug, lattice, IO, crystallite ! diff --git a/trunk/mpie_cpfem_marc2007r1_sequential.f90 b/trunk/mpie_cpfem_marc2007r1_sequential.f90 index 313de55df..5a090d833 100644 --- a/trunk/mpie_cpfem_marc2007r1_sequential.f90 +++ b/trunk/mpie_cpfem_marc2007r1_sequential.f90 @@ -4,7 +4,7 @@ ! written by F. Roters, P. Eisenlohr, L. Hantcherli, W.A. Counts ! MPI fuer Eisenforschung, Duesseldorf ! -! last modified: 09.07.2008 +! last modified: 27.11.2008 !******************************************************************** ! Usage: ! - choose material as hypela2 @@ -35,7 +35,7 @@ include "mesh.f90" ! uses prec, IO, math, FEsolving include "lattice.f90" ! uses prec, math include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug - include "crystallite.f90" ! uses prec, debug, constitutive, mesh, math, IO +! include "crystallite.f90" ! uses prec, debug, constitutive, mesh, math, IO include "CPFEM_sequential.f90" ! uses prec, math, mesh, constitutive, FEsolving, debug, lattice, IO, crystallite ! diff --git a/trunk/mpie_cpfem_marc2008r1.f90 b/trunk/mpie_cpfem_marc2008r1.f90 index d042fcc57..d66e407cd 100644 --- a/trunk/mpie_cpfem_marc2008r1.f90 +++ b/trunk/mpie_cpfem_marc2008r1.f90 @@ -4,7 +4,7 @@ ! written by F. Roters, P. Eisenlohr, L. Hantcherli, W.A. Counts ! MPI fuer Eisenforschung, Duesseldorf ! -! last modified: 28.10.2008 +! last modified: 22.11.2008 !******************************************************************** ! Usage: ! - choose material as hypela2 @@ -36,7 +36,7 @@ include "mesh.f90" ! uses prec, IO, math, FEsolving include "lattice.f90" ! uses prec, math include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug - include "crystallite.f90" ! uses prec, debug, constitutive, mesh, math, IO +! include "crystallite.f90" ! uses prec, debug, constitutive, mesh, math, IO include "CPFEM.f90" ! uses prec, math, mesh, constitutive, FEsolving, debug, lattice, IO, crystallite ! ! diff --git a/trunk/mpie_cpfem_marc2008r1_sequential.f90 b/trunk/mpie_cpfem_marc2008r1_sequential.f90 index 155da786c..0f9e165ee 100644 --- a/trunk/mpie_cpfem_marc2008r1_sequential.f90 +++ b/trunk/mpie_cpfem_marc2008r1_sequential.f90 @@ -4,7 +4,7 @@ ! written by F. Roters, P. Eisenlohr, L. Hantcherli, W.A. Counts ! MPI fuer Eisenforschung, Duesseldorf ! -! last modified: 09.07.2008 +! last modified: 27.11.2008 !******************************************************************** ! Usage: ! - choose material as hypela2 @@ -16,6 +16,7 @@ ! - use nonsymmetric option for solver (e.g. direct ! profile or multifrontal sparse, the latter seems ! to be faster!) +! - in case of ddm a symmetric solver has to be used !******************************************************************** ! Marc subroutines used: ! - hypela2 @@ -35,7 +36,7 @@ include "mesh.f90" ! uses prec, IO, math, FEsolving include "lattice.f90" ! uses prec, math include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug - include "crystallite.f90" ! uses prec, debug, constitutive, mesh, math, IO +! include "crystallite.f90" ! uses prec, debug, constitutive, mesh, math, IO include "CPFEM_sequential.f90" ! uses prec, math, mesh, constitutive, FEsolving, debug, lattice, IO, crystallite ! diff --git a/trunk/todo b/trunk/todo.txt similarity index 73% rename from trunk/todo rename to trunk/todo.txt index b412619f8..b4ffd9d03 100644 --- a/trunk/todo +++ b/trunk/todo.txt @@ -1,5 +1,8 @@ -Things to be implemented into the code - -# parsing of texture data from mattex file to be done by separate "texture.f90", thus freeing constitutive.f90 from this global task -# constitutive law to be used is indicated within [material]. constitutive.f90 then assumes a wrapper functionality for the four functions "dor_microstructure", "Lpanditstangent", "postresults" and "microstructure". Within those, a switch case checks for the true constitutive law (a string) and passes control on to the respective subroutine-clone of this law. Technically, the compilation relies on "IFDEF" to include the necessary constitutive laws plus the corresponding entries in each switch case. -# change state variable meaning to (i) homogenization, (ii) microstructure +Things to be implemented into the code + +# parsing of texture data from mattex file to be done by separate "texture.f90", thus freeing constitutive.f90 from this global task + this should go along with new scheme for multi material definition +# constitutive law to be used is indicated within [material]. constitutive.f90 then assumes a wrapper functionality for the four functions "dot_microstructure", "Lpanditstangent", "postresults" and "microstructure". Within those, a switch case checks for the true constitutive law (a string) and passes control on to the respective subroutine-clone of this law. Technically, the compilation relies on "IFDEF" to include the necessary constitutive laws plus the corresponding entries in each switch case. +# change state variable meaning to (i) homogenization, (ii) microstructure +# adopt CPFEM_GIA8.f90 to new scheme, rename to CPFEM_RGC.f90 +# make OpenMP parallelization work again