diff --git a/CMakeLists.txt b/CMakeLists.txt index a5e021a2b..b8d344240 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -116,7 +116,7 @@ set (DAMASK_VERSION_MINOR ${DAMASK_V}) # Built-in options for DAMASK build system # -> can be overwritten from commandline/install_script option(OPENMP "Use OpenMP libaries for DAMASK" ON ) -option(OPTIMIZATION "DAMASK optimization level [OFF,DEFENSIVE,AGGRESSIVE]" "DEFENSIVE" ) +option(OPTIMIZATION "DAMASK optimization level [OFF,DEFENSIVE,AGGRESSIVE]" "AGGRESSIVE" ) option(SPECTRAL "Build spectral sovler for DAMASAK" OFF ) option(FEM "Build FEM solver for DAMASK" OFF ) diff --git a/Makefile b/Makefile index 5c4caa4a5..24b17965f 100755 --- a/Makefile +++ b/Makefile @@ -10,7 +10,7 @@ spectral: build/spectral build/spectral: build @mkdir build/spectral - @(cd build/spectral; cmake -Wno-dev -DCMAKE_BUILD_TYPE=RELEASE -DDAMASK_DRIVER=SPECTRAL ../..;) + @(cd build/spectral; cmake -Wno-dev -DCMAKE_BUILD_TYPE=RELEASE -DDAMASK_DRIVER=SPECTRAL -DOPTIMIZATION=AGGRESSIVE -DDAMASK_INSTALL=${HOME}/bin ../..;) build: bin @mkdir build diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 25e56f4ed..9e5875629 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -79,7 +79,6 @@ add_library(DAMASK_DRIVERS ALIAS DAMASK_LATTICE) add_library (DAMASK_PLASTIC "plastic_dislotwin.f90" "plastic_disloUCLA.f90" "plastic_isotropic.f90" - "plastic_j2.f90" "plastic_phenopowerlaw.f90" "plastic_titanmod.f90" "plastic_nonlocal.f90" diff --git a/src/constitutive.f90 b/src/constitutive.f90 index c72af06ca..3ef51aecf 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -207,7 +207,7 @@ subroutine constitutive_init() outputName = PLASTICITY_NONE_label thisNoutput => null() thisOutput => null() - thisSize => null() + thisSize => null() case (PLASTICITY_ISOTROPIC_ID) plasticityType outputName = PLASTICITY_ISOTROPIC_label thisNoutput => plastic_isotropic_Noutput @@ -423,7 +423,7 @@ end function constitutive_homogenizedC !-------------------------------------------------------------------------------------------------- !> @brief calls microstructure function of the different constitutive models !-------------------------------------------------------------------------------------------------- -subroutine constitutive_microstructure(orientations, Fe, Fp, ipc, ip, el) +subroutine constitutive_microstructure(orientations, Fe, Fp, ipc, ip, el, F0s,Fes,Fps,Tstar_vs) use prec, only: & pReal use material, only: & @@ -460,7 +460,15 @@ subroutine constitutive_microstructure(orientations, Fe, Fp, ipc, ip, el) ho, & !< homogenization tme !< thermal member position real(pReal), intent(in), dimension(:,:,:,:) :: & - orientations !< crystal orientations as quaternions + orientations + + real(pReal), intent(in), dimension(:,:,:,:,:) :: & + F0s, & + Fes, & + Fps + + real(pReal), intent(in), dimension(:,:,:,:) :: & + Tstar_vs !< crystal orientations as quaternions ho = material_homog(ip,el) tme = thermalMapping(ho)%p(ip,el) @@ -475,7 +483,7 @@ subroutine constitutive_microstructure(orientations, Fe, Fp, ipc, ip, el) case (PLASTICITY_NONLOCAL_ID) plasticityType call plastic_nonlocal_microstructure (Fe,Fp,ip,el) case (PLASTICITY_PHENOPLUS_ID) plasticityType - call plastic_phenoplus_microstructure(orientations,ipc,ip,el) + call plastic_phenoplus_microstructure(orientations,ipc,ip,el,F0s,Fes,Fps,Tstar_vs) end select plasticityType end subroutine constitutive_microstructure diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 60fc373b7..98159d7a3 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -425,7 +425,7 @@ subroutine crystallite_init crystallite_partionedFp0 = crystallite_Fp0 crystallite_partionedFi0 = crystallite_Fi0 crystallite_partionedF0 = crystallite_F0 - crystallite_partionedF = crystallite_F0 + crystallite_partionedF = crystallite_F0 call crystallite_orientations() crystallite_orientation0 = crystallite_orientation ! store initial orientations for calculation of grain rotations @@ -438,7 +438,11 @@ subroutine crystallite_init call constitutive_microstructure(crystallite_orientation, & ! pass orientation to constitutive module crystallite_Fe(1:3,1:3,c,i,e), & crystallite_Fp(1:3,1:3,c,i,e), & - c,i,e) ! update dependent state variables to be consistent with basic states + c,i,e, & + crystallite_F0, & + crystallite_Fe, & + crystallite_Fp, & + crystallite_Tstar_v) ! update dependent state variables to be consistent with basic states enddo enddo enddo @@ -1715,7 +1719,11 @@ subroutine crystallite_integrateStateRK4() call constitutive_microstructure(crystallite_orientation, & crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), & - g, i, e) ! update dependent state variables to be consistent with basic states + g, i, e, & + crystallite_F0, & + crystallite_Fe, & + crystallite_Fp, & + crystallite_Tstar_v) ! update dependent state variables to be consistent with basic states enddo; enddo; enddo !$OMP ENDDO @@ -2041,7 +2049,11 @@ subroutine crystallite_integrateStateRKCK45() call constitutive_microstructure(crystallite_orientation, & crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), & - g, i, e) ! update dependent state variables to be consistent with basic states + g, i, e, & + crystallite_F0, & + crystallite_Fe, & + crystallite_Fp, & + crystallite_Tstar_v) ! update dependent state variables to be consistent with basic states enddo; enddo; enddo !$OMP ENDDO @@ -2261,7 +2273,11 @@ subroutine crystallite_integrateStateRKCK45() call constitutive_microstructure(crystallite_orientation, & crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), & - g, i, e) ! update dependent state variables to be consistent with basic states + g, i, e, & + crystallite_F0, & + crystallite_Fe, & + crystallite_Fp, & + crystallite_Tstar_v) ! update dependent state variables to be consistent with basic states enddo; enddo; enddo !$OMP ENDDO @@ -2496,7 +2512,11 @@ subroutine crystallite_integrateStateAdaptiveEuler() call constitutive_microstructure(crystallite_orientation, & crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), & - g, i, e) ! update dependent state variables to be consistent with basic states + g, i, e, & + crystallite_F0, & + crystallite_Fe, & + crystallite_Fp, & + crystallite_Tstar_v) ! update dependent state variables to be consistent with basic states enddo; enddo; enddo !$OMP ENDDO !$OMP END PARALLEL @@ -2840,7 +2860,11 @@ eIter = FEsolving_execElem(1:2) call constitutive_microstructure(crystallite_orientation, & crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), & - g, i, e) ! update dependent state variables to be consistent with basic states + g, i, e, & + crystallite_F0, & + crystallite_Fe, & + crystallite_Fp, & + crystallite_Tstar_v) ! update dependent state variables to be consistent with basic states enddo; enddo; enddo !$OMP ENDDO !$OMP END PARALLEL @@ -3085,7 +3109,11 @@ subroutine crystallite_integrateStateFPI() call constitutive_microstructure(crystallite_orientation, & crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), & - g, i, e) ! update dependent state variables to be consistent with basic states + g, i, e, & + crystallite_F0, & + crystallite_Fe, & + crystallite_Fp, & + crystallite_Tstar_v) ! update dependent state variables to be consistent with basic states p = phaseAt(g,i,e) c = phasememberAt(g,i,e) plasticState(p)%previousDotState2(:,c) = plasticState(p)%previousDotState(:,c) diff --git a/src/lattice.f90 b/src/lattice.f90 index bf1f53658..4ba915688 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -18,7 +18,7 @@ module lattice LATTICE_maxNtwinFamily = 4_pInt, & !< max # of twin system families over lattice structures LATTICE_maxNtransFamily = 2_pInt, & !< max # of transformation system families over lattice structures LATTICE_maxNcleavageFamily = 3_pInt !< max # of transformation system families over lattice structures - + integer(pInt), allocatable, dimension(:,:), protected, public :: & lattice_NslipSystem, & !< total # of slip systems in each family lattice_NtwinSystem, & !< total # of twin systems in each family @@ -26,22 +26,22 @@ module lattice lattice_NcleavageSystem !< total # of transformation systems in each family integer(pInt), allocatable, dimension(:,:,:), protected, public :: & - lattice_interactionSlipSlip, & !< Slip--slip interaction type - lattice_interactionSlipTwin, & !< Slip--twin interaction type - lattice_interactionTwinSlip, & !< Twin--slip interaction type - lattice_interactionTwinTwin, & !< Twin--twin interaction type - lattice_interactionSlipTrans, & !< Slip--trans interaction type - lattice_interactionTransSlip, & !< Trans--slip interaction type + lattice_interactionSlipSlip, & !< Slip--slip interaction type + lattice_interactionSlipTwin, & !< Slip--twin interaction type + lattice_interactionTwinSlip, & !< Twin--slip interaction type + lattice_interactionTwinTwin, & !< Twin--twin interaction type + lattice_interactionSlipTrans, & !< Slip--trans interaction type + lattice_interactionTransSlip, & !< Trans--slip interaction type lattice_interactionTransTrans !< Trans--trans interaction type real(pReal), allocatable, dimension(:,:,:,:,:), protected, public :: & lattice_Sslip, & !< Schmid and non-Schmid matrices lattice_Scleavage !< Schmid matrices for cleavage systems - + real(pReal), allocatable, dimension(:,:,:,:), protected, public :: & lattice_Sslip_v, & !< Mandel notation of lattice_Sslip lattice_Scleavage_v !< Mandel notation of lattice_Scleavege - + real(pReal), allocatable, dimension(:,:,:), protected, public :: & lattice_sn, & !< normal direction of slip system lattice_sd, & !< slip direction of slip system @@ -75,25 +75,25 @@ module lattice !-------------------------------------------------------------------------------------------------- ! face centered cubic - integer(pInt), dimension(LATTICE_maxNslipFamily), parameter, public :: & + integer(pInt), dimension(LATTICE_maxNslipFamily), parameter, public :: & LATTICE_fcc_NslipSystem = int([12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],pInt) !< # of slip systems per family for fcc - + integer(pInt), dimension(LATTICE_maxNtwinFamily), parameter, public :: & LATTICE_fcc_NtwinSystem = int([12, 0, 0, 0],pInt) !< # of twin systems per family for fcc integer(pInt), dimension(LATTICE_maxNtransFamily), parameter, public :: & LATTICE_fcc_NtransSystem = int([12, 0],pInt) !< # of transformation systems per family for fcc - + integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: & LATTICE_fcc_NcleavageSystem = int([3, 4, 0],pInt) !< # of cleavage systems per family for fcc - + integer(pInt), parameter, private :: & LATTICE_fcc_Nslip = 12_pInt, & !sum(lattice_fcc_NslipSystem), & !< total # of slip systems for fcc LATTICE_fcc_Ntwin = 12_pInt, & !sum(lattice_fcc_NtwinSystem), & !< total # of twin systems for fcc LATTICE_fcc_NnonSchmid = 0_pInt, & !< total # of non-Schmid contributions for fcc LATTICE_fcc_Ntrans = 12_pInt, & !sum(lattice_fcc_NtransSystem), & !< total # of transformation systems for fcc LATTICE_fcc_Ncleavage = 7_pInt !sum(lattice_fcc_NcleavageSystem) !< total # of cleavage systems for fcc - + real(pReal), dimension(3+3,LATTICE_fcc_Nslip), parameter, private :: & LATTICE_fcc_systemSlip = reshape(real([& ! Slip direction Plane normal @@ -254,10 +254,10 @@ module lattice 2,2,2,2,2,2,2,2,2,1,1,1, & 2,2,2,2,2,2,2,2,2,1,1,1 & ],pInt),[LATTICE_fcc_Ntrans,LATTICE_fcc_Ntrans],order=[2,1]) !< Trans--trans interaction types for fcc - + real(pReal), dimension(LATTICE_fcc_Ntrans), parameter, private :: & LATTICE_fccTohex_shearTrans = sqrt(2.0_pReal)/4.0_pReal - + real(pReal), dimension(4,LATTICE_fcc_Ntrans), parameter, private :: & LATTICE_fccTobcc_systemTrans = reshape([& 0.0, 1.0, 0.0, 10.26, & ! Pitsch OR (Ma & Hartmaier 2014, Table 3) @@ -287,7 +287,7 @@ module lattice 0, 0, 1, 1, 0, 0, 0, 1, 0, & 0, 0, 1, 1, 0, 0, 0, 1, 0, & 0, 0, 1, 1, 0, 0, 0, 1, 0, & - 0, 0, 1, 1, 0, 0, 0, 1, 0 & + 0, 0, 1, 1, 0, 0, 0, 1, 0 & ],pInt),[ 9_pInt, LATTICE_fcc_Ntrans]) real(pReal), dimension(4,LATTICE_fcc_Ntrans), parameter, private :: & @@ -308,7 +308,7 @@ module lattice real(pReal), dimension(LATTICE_fcc_Ntrans,LATTICE_fcc_Ntrans), parameter, private :: & ! Matrix for projection of shear from slip system to fault-band (twin) systems LATTICE_fccTobcc_projectionTrans = reshape(real([& ! For ns = nt = nr - 0, 1,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 0, 1,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, & -1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, & 1,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & 0, 0, 0, 0, 1,-1, 0, 0, 0, 0, 0, 0, & @@ -357,10 +357,10 @@ module lattice ],pReal),[ 3_pInt + 3_pInt,LATTICE_fcc_Ncleavage]) !-------------------------------------------------------------------------------------------------- -! body centered cubic +! body centered cubic integer(pInt), dimension(LATTICE_maxNslipFamily), parameter, public :: & LATTICE_bcc_NslipSystem = int([ 12, 12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], pInt) !< # of slip systems per family for bcc - + integer(pInt), dimension(LATTICE_maxNtwinFamily), parameter, public :: & LATTICE_bcc_NtwinSystem = int([ 12, 0, 0, 0], pInt) !< # of twin systems per family for bcc @@ -369,7 +369,7 @@ module lattice integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: & LATTICE_bcc_NcleavageSystem = int([3,6,0],pInt) !< # of cleavage systems per family for bcc - + integer(pInt), parameter, private :: & LATTICE_bcc_Nslip = 24_pInt, & !sum(lattice_bcc_NslipSystem), & !< total # of slip systems for bcc LATTICE_bcc_Ntwin = 12_pInt, & !sum(lattice_bcc_NtwinSystem), & !< total # of twin systems for bcc @@ -377,10 +377,11 @@ module lattice LATTICE_bcc_Ntrans = 0_pInt, & !sum(lattice_bcc_NtransSystem), & !< total # of transformation systems for bcc LATTICE_bcc_Ncleavage = 9_pInt !sum(lattice_bcc_NcleavageSystem) !< total # of cleavage systems for bcc + real(pReal), dimension(3+3,LATTICE_bcc_Nslip), parameter, private :: & LATTICE_bcc_systemSlip = reshape(real([& ! Slip direction Plane normal - ! Slip system <111>{110} + ! Slip system <111>{110} 1,-1, 1, 0, 1, 1, & -1,-1, 1, 0, 1, 1, & 1, 1, 1, 0,-1, 1, & @@ -455,7 +456,7 @@ module lattice integer(pInt), dimension(LATTICE_bcc_Nslip,LATTICE_bcc_Nslip), parameter, public :: & LATTICE_bcc_interactionSlipSlip = reshape(int( [& - 1,2,6,6,5,4,4,3,4,3,5,4, 6,6,4,3,3,4,6,6,4,3,6,6, & ! ---> slip + 1,2,6,6,5,4,4,3,4,3,5,4, 6,6,4,3,3,4,6,6,4,3,6,6, & ! ---> slip 2,1,6,6,4,3,5,4,5,4,4,3, 6,6,3,4,4,3,6,6,3,4,6,6, & ! | 6,6,1,2,4,5,3,4,4,5,3,4, 4,3,6,6,6,6,3,4,6,6,4,3, & ! | 6,6,2,1,3,4,4,5,3,4,4,5, 3,4,6,6,6,6,4,3,6,6,3,4, & ! v slip @@ -492,7 +493,7 @@ module lattice 3,3,3,2,2,3,3,3,3,2,3,3, & ! ---> twin 3,3,2,3,3,2,3,3,2,3,3,3, & ! | 3,2,3,3,3,3,2,3,3,3,3,2, & ! | - 2,3,3,3,3,3,3,2,3,3,2,3, & ! v slip + 2,3,3,3,3,3,3,2,3,3,2,3, & ! v slip 2,3,3,3,3,3,3,2,3,3,2,3, & 3,3,2,3,3,2,3,3,2,3,3,3, & 3,2,3,3,3,3,2,3,3,3,3,2, & @@ -556,17 +557,17 @@ module lattice !-------------------------------------------------------------------------------------------------- ! hexagonal integer(pInt), dimension(LATTICE_maxNslipFamily), parameter, public :: & - lattice_hex_NslipSystem = int([ 3, 3, 3, 6, 12, 6, 0, 0, 0, 0, 0, 0, 0],pInt) !< # of slip systems per family for hex - + lattice_hex_NslipSystem = int([ 3, 3, 3, 6, 12, 6, 0, 0, 0, 0, 0, 0, 0],pInt) !< # of slip systems per family for hex + integer(pInt), dimension(LATTICE_maxNtwinFamily), parameter, public :: & - lattice_hex_NtwinSystem = int([ 6, 6, 6, 6],pInt) !< # of slip systems per family for hex + lattice_hex_NtwinSystem = int([ 6, 6, 6, 6],pInt) !< # of slip systems per family for hex integer(pInt), dimension(LATTICE_maxNtransFamily), parameter, public :: & LATTICE_hex_NtransSystem = int([0,0],pInt) !< # of transformation systems per family for hex - + integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: & LATTICE_hex_NcleavageSystem = int([3,0,0],pInt) !< # of cleavage systems per family for hex - + integer(pInt), parameter, private :: & LATTICE_hex_Nslip = 33_pInt, & !sum(lattice_hex_NslipSystem), & !< total # of slip systems for hex LATTICE_hex_Ntwin = 24_pInt, & !sum(lattice_hex_NtwinSystem), & !< total # of twin systems for hex @@ -609,7 +610,7 @@ module lattice 2, -1, -1, 3, -1, 0, 1, 1, & 1, -2, 1, 3, 0, 1, -1, 1, & -1, -1, 2, 3, 0, 1, -1, 1, & - ! pyramidal system: c+a slip <11.3>{-1-1.2} -- as for hexagonal ice (Castelnau et al. 1996, similar to twin system found below) + ! pyramidal system: c+a slip <11.3>{-1-1.2} -- as for hexagonal ice (Castelnau et al. 1996, similar to twin system found below) 2, -1, -1, 3, -2, 1, 1, 2, & ! sorted according to similar twin system -1, 2, -1, 3, 1, -2, 1, 2, & ! <11.3>{-1-1.2} shear = 2((c/a)^2-2)/(3 c/a) -1, -1, 2, 3, 1, 1, -2, 2, & @@ -648,7 +649,7 @@ module lattice -2, 1, 1, -3, -2, 1, 1, 2, & 1, -2, 1, -3, 1, -2, 1, 2, & 1, 1, -2, -3, 1, 1, -2, 2 & - ],pReal),[ 4_pInt + 4_pInt ,LATTICE_hex_Ntwin]) !< twin systems for hex, order follows Prof. Tom Bieler's scheme; but numbering in data was restarted from 1 + ],pReal),[ 4_pInt + 4_pInt ,LATTICE_hex_Ntwin]) !< twin systems for hex, order follows Prof. Tom Bieler's scheme; but numbering in data was restarted from 1 integer(pInt), dimension(LATTICE_hex_Ntwin), parameter, private :: & LATTICE_hex_shearTwin = reshape(int( [& ! indicator to formula further below @@ -677,7 +678,7 @@ module lattice 4, & 4 & ],pInt),[LATTICE_hex_Ntwin]) - + integer(pInt), dimension(LATTICE_hex_Nslip,LATTICE_hex_Nslip), parameter, public :: & LATTICE_hex_interactionSlipSlip = reshape(int( [& 1, 2, 2, 3, 3, 3, 7, 7, 7, 13,13,13,13,13,13, 21,21,21,21,21,21,21,21,21,21,21,21, 31,31,31,31,31,31, & ! ---> slip @@ -717,10 +718,10 @@ module lattice 42,42,42, 41,41,41, 40,40,40, 39,39,39,39,39,39, 38,38,38,38,38,38,38,38,38,38,38,38, 37,37,36,37,37,37, & 42,42,42, 41,41,41, 40,40,40, 39,39,39,39,39,39, 38,38,38,38,38,38,38,38,38,38,38,38, 37,37,37,36,37,37, & 42,42,42, 41,41,41, 40,40,40, 39,39,39,39,39,39, 38,38,38,38,38,38,38,38,38,38,38,38, 37,37,37,37,36,37, & - 42,42,42, 41,41,41, 40,40,40, 39,39,39,39,39,39, 38,38,38,38,38,38,38,38,38,38,38,38, 37,37,37,37,37,36 & + 42,42,42, 41,41,41, 40,40,40, 39,39,39,39,39,39, 38,38,38,38,38,38,38,38,38,38,38,38, 37,37,37,37,37,36 & ! ],pInt),[LATTICE_hex_Nslip,LATTICE_hex_Nslip],order=[2,1]) !< Slip--slip interaction types for hex (onion peel naming scheme) - + integer(pInt), dimension(LATTICE_hex_Nslip,LATTICE_hex_Ntwin), parameter, public :: & LATTICE_hex_interactionSlipTwin = reshape(int( [& 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, & ! --> twin @@ -728,7 +729,7 @@ module lattice 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, & ! | ! v 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, & ! slip - 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, & + 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, & 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, & ! 9, 9, 9, 9, 9, 9, 10,10,10,10,10,10, 11,11,11,11,11,11, 12,12,12,12,12,12, & @@ -762,7 +763,7 @@ module lattice 21,21,21,21,21,21, 22,22,22,22,22,22, 23,23,23,23,23,23, 24,24,24,24,24,24, & 21,21,21,21,21,21, 22,22,22,22,22,22, 23,23,23,23,23,23, 24,24,24,24,24,24 & ! - ],pInt),[LATTICE_hex_Nslip,LATTICE_hex_Ntwin],order=[2,1]) !< Slip--twin interaction types for hex (isotropic, 24 in total) + ],pInt),[LATTICE_hex_Nslip,LATTICE_hex_Ntwin],order=[2,1]) !< Slip--twin interaction types for hex (isotropic, 24 in total) integer(pInt), dimension(LATTICE_hex_Ntwin,LATTICE_hex_Nslip), parameter, public :: & LATTICE_hex_interactionTwinSlip = reshape(int( [& @@ -835,6 +836,7 @@ module lattice ],pReal),[ 4_pInt + 4_pInt,LATTICE_hex_Ncleavage]) + !-------------------------------------------------------------------------------------------------- ! body centered tetragonal integer(pInt), dimension(LATTICE_maxNslipFamily), parameter, public :: & @@ -848,14 +850,14 @@ module lattice integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: & LATTICE_bct_NcleavageSystem = int([0,0,0],pInt) !< # of cleavage systems per family for bct - + integer(pInt), parameter, private :: & LATTICE_bct_Nslip = 52_pInt, & !sum(lattice_bct_NslipSystem), & !< total # of slip systems for bct LATTICE_bct_Ntwin = 0_pInt, & !sum(lattice_bct_NtwinSystem), & !< total # of twin systems for bct LATTICE_bct_NnonSchmid = 0_pInt, & !< total # of non-Schmid contributions for bct LATTICE_bct_Ntrans = 0_pInt, & !sum(lattice_bct_NtransSystem), & !< total # of transformation systems for bct LATTICE_bct_Ncleavage = 0_pInt !sum(lattice_bct_NcleavageSystem) !< total # of cleavage systems for bct - + real(pReal), dimension(3+3,LATTICE_bct_Nslip), parameter, private :: & LATTICE_bct_systemSlip = reshape(real([& ! Slip direction Plane normal @@ -865,7 +867,7 @@ module lattice ! Slip family 2 {110)<001] 0, 0, 1, 1, 1, 0, & 0, 0, 1, -1, 1, 0, & - ! slip family 3 {100)<010] + ! slip family 3 {100)<010] 0, 1, 0, 1, 0, 0, & 1, 0, 0, 0, 1, 0, & ! Slip family 4 {110)<1-11]/2 @@ -876,54 +878,54 @@ module lattice ! Slip family 5 {110)<1-10] 1, -1, 0, 1, 1, 0, & 1, 1, 0, 1,-1, 0, & - ! Slip family 6 {100)<011] + ! Slip family 6 {100)<011] 0, 1, 1, 1, 0, 0, & 0,-1, 1, 1, 0, 0, & -1, 0, 1, 0, 1, 0, & - 1, 0, 1, 0, 1, 0, & - ! Slip family 7 {001)<010] + 1, 0, 1, 0, 1, 0, & + ! Slip family 7 {001)<010] 0, 1, 0, 0, 0, 1, & 1, 0, 0, 0, 0, 1, & - ! Slip family 8 {001)<110] + ! Slip family 8 {001)<110] 1, 1, 0, 0, 0, 1, & - -1, 1, 0, 0, 0, 1, & - ! Slip family 9 {011)<01-1] + -1, 1, 0, 0, 0, 1, & + ! Slip family 9 {011)<01-1] 0, 1,-1, 0, 1, 1, & 0,-1,-1, 0,-1, 1, & -1, 0,-1, -1, 0, 1, & - 1, 0,-1, 1, 0, 1, & - ! Slip family 10 {011)<1-11]/2 + 1, 0,-1, 1, 0, 1, & + ! Slip family 10 {011)<1-11]/2 1,-1, 1, 0, 1, 1, & 1, 1,-1, 0, 1, 1, & 1, 1, 1, 0, 1,-1, & - -1, 1, 1, 0, 1,-1, & + -1, 1, 1, 0, 1,-1, & 1,-1,-1, 1, 0, 1, & -1,-1, 1, 1, 0, 1, & 1, 1, 1, 1, 0,-1, & - 1,-1, 1, 1, 0,-1, & - ! Slip family 11 {011)<100] + 1,-1, 1, 1, 0,-1, & + ! Slip family 11 {011)<100] 1, 0, 0, 0, 1, 1, & 1, 0, 0, 0, 1,-1, & 0, 1, 0, 1, 0, 1, & - 0, 1, 0, 1, 0,-1, & - ! Slip family 12 {211)<01-1] + 0, 1, 0, 1, 0,-1, & + ! Slip family 12 {211)<01-1] 0, 1,-1, 2, 1, 1, & 0,-1,-1, 2,-1, 1, & 1, 0,-1, 1, 2, 1, & - -1, 0,-1, -1, 2, 1, & + -1, 0,-1, -1, 2, 1, & 0, 1,-1, -2, 1, 1, & 0,-1,-1, -2,-1, 1, & -1, 0,-1, -1,-2, 1, & - 1, 0,-1, 1,-2, 1, & - ! Slip family 13 {211)<-111]/2 + 1, 0,-1, 1,-2, 1, & + ! Slip family 13 {211)<-111]/2 -1, 1, 1, 2, 1, 1, & -1,-1, 1, 2,-1, 1, & 1,-1, 1, 1, 2, 1, & - -1,-1, 1, -1, 2, 1, & + -1,-1, 1, -1, 2, 1, & 1, 1, 1, -2, 1, 1, & 1,-1, 1, -2,-1, 1, & -1, 1, 1, -1,-2, 1, & - 1, 1, 1, 1,-2, 1 & + 1, 1, 1, 1,-2, 1 & ],pReal),[ 3_pInt + 3_pInt,LATTICE_bct_Nslip]) !< slip systems for bct sorted by Bieler integer(pInt), dimension(LATTICE_bct_Nslip,LATTICE_bct_Nslip), parameter, public :: & @@ -933,7 +935,7 @@ module lattice ! 6, 6, 4, 5, 8, 8, 14, 14, 14, 14, 22, 22, 32, 32, 32, 32, 44, 44, 58, 58, 74, 74, 74, 74, 92, 92, 92, 92, 92, 92, 92, 92, 112, 112, 112, 112, 134,134,134,134,134,134,134,134, 158,158,158,158,158,158,158,158, & 6, 6, 5, 4, 8, 8, 14, 14, 14, 14, 22, 22, 32, 32, 32, 32, 44, 44, 58, 58, 74, 74, 74, 74, 92, 92, 92, 92, 92, 92, 92, 92, 112, 112, 112, 112, 134,134,134,134,134,134,134,134, 158,158,158,158,158,158,158,158, & - ! + ! 12, 12, 11, 11, 9, 10, 15, 15, 15, 15, 23, 23, 33, 33, 33, 33, 45, 45, 59, 59, 75, 75, 75, 75, 93, 93, 93, 93, 93, 93, 93, 93, 113, 113, 113, 113, 135,135,135,135,135,135,135,135, 159,159,159,159,159,159,159,159, & 12, 12, 11, 11, 10, 9, 15, 15, 15, 15, 23, 23, 33, 33, 33, 33, 45, 45, 59, 59, 75, 75, 75, 75, 93, 93, 93, 93, 93, 93, 93, 93, 113, 113, 113, 113, 135,135,135,135,135,135,135,135, 159,159,159,159,159,159,159,159, & ! @@ -993,7 +995,7 @@ module lattice 182,182, 181,181, 180,180, 179,179,179,179, 178,178, 177,177,177,177, 176,176, 175,175, 174,174,174,174, 173,173,173,173,173,173,173,173, 172, 172, 172, 172, 171,171,171,171,171,171,171,171, 169,170,170,170,170,170,169,170, & 182,182, 181,181, 180,180, 179,179,179,179, 178,178, 177,177,177,177, 176,176, 175,175, 174,174,174,174, 173,173,173,173,173,173,173,173, 172, 172, 172, 172, 171,171,171,171,171,171,171,171, 169,170,170,170,170,170,170,169 & - ],pInt),[lattice_bct_Nslip,lattice_bct_Nslip],order=[2,1]) + ],pInt),[lattice_bct_Nslip,lattice_bct_Nslip],order=[2,1]) !-------------------------------------------------------------------------------------------------- ! isotropic @@ -1008,14 +1010,14 @@ module lattice integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: & LATTICE_iso_NcleavageSystem = int([3,0,0],pInt) !< # of cleavage systems per family for iso - + integer(pInt), parameter, private :: & LATTICE_iso_Nslip = 0_pInt, & !sum(lattice_iso_NslipSystem), & !< total # of slip systems for iso LATTICE_iso_Ntwin = 0_pInt, & !sum(lattice_iso_NtwinSystem), & !< total # of twin systems for iso LATTICE_iso_NnonSchmid = 0_pInt, & !< total # of non-Schmid contributions for iso LATTICE_iso_Ntrans = 0_pInt, & !sum(lattice_iso_NtransSystem), & !< total # of transformation systems for iso LATTICE_iso_Ncleavage = 3_pInt !sum(lattice_iso_NcleavageSystem) !< total # of cleavage systems for iso - + real(pReal), dimension(3+3,LATTICE_iso_Ncleavage), parameter, private :: & LATTICE_iso_systemCleavage = reshape(real([& ! Cleavage direction Plane normal @@ -1037,14 +1039,14 @@ module lattice integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: & LATTICE_ortho_NcleavageSystem = int([1,1,1],pInt) !< # of cleavage systems per family for ortho - + integer(pInt), parameter, private :: & LATTICE_ortho_Nslip = 0_pInt, & !sum(lattice_ortho_NslipSystem), & !< total # of slip systems for ortho LATTICE_ortho_Ntwin = 0_pInt, & !sum(lattice_ortho_NtwinSystem), & !< total # of twin systems for ortho LATTICE_ortho_NnonSchmid = 0_pInt, & !< total # of non-Schmid contributions for ortho LATTICE_ortho_Ntrans = 0_pInt, & !sum(lattice_ortho_NtransSystem), & !< total # of transformation systems for ortho LATTICE_ortho_Ncleavage = 3_pInt !sum(lattice_ortho_NcleavageSystem) !< total # of cleavage systems for ortho - + real(pReal), dimension(3+3,LATTICE_ortho_Ncleavage), parameter, private :: & LATTICE_ortho_systemCleavage = reshape(real([& ! Cleavage direction Plane normal @@ -1075,14 +1077,14 @@ module lattice real(pReal), dimension(:,:,:), allocatable, public, protected :: & lattice_C66, lattice_trans_C66 - real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: & + real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: & lattice_C3333, lattice_trans_C3333 - real(pReal), dimension(:), allocatable, public, protected :: & + real(pReal), dimension(:), allocatable, public, protected :: & lattice_mu, & lattice_nu, & lattice_trans_mu, & lattice_trans_nu - real(pReal), dimension(:,:,:), allocatable, public, protected :: & + real(pReal), dimension(:,:,:), allocatable, public, protected :: & lattice_thermalConductivity33, & lattice_thermalExpansion33, & lattice_damageDiffusion33, & @@ -1091,7 +1093,7 @@ module lattice lattice_porosityDiffusion33, & lattice_hydrogenfluxDiffusion33, & lattice_hydrogenfluxMobility33 - real(pReal), dimension(:), allocatable, public, protected :: & + real(pReal), dimension(:), allocatable, public, protected :: & lattice_damageMobility, & lattice_porosityMobility, & lattice_massDensity, & @@ -1120,7 +1122,7 @@ module lattice integer(pInt), dimension(2), parameter, private :: & - lattice_NsymOperations = [24_pInt,12_pInt] + lattice_NsymOperations = [24_pInt,12_pInt] real(pReal), dimension(4,36), parameter, private :: & lattice_symOperations = reshape([& @@ -1268,7 +1270,7 @@ subroutine lattice_init debug_levelBasic use numerics, only: & worldrank - + implicit none integer(pInt), parameter :: FILEUNIT = 200_pInt integer(pInt) :: Nphases @@ -1279,11 +1281,11 @@ subroutine lattice_init integer(pInt) :: section = 0_pInt,i real(pReal), dimension(:), allocatable :: & CoverA, & !!!!!!< c/a ratio for low symmetry type lattice - CoverA_trans, & !< c/a ratio for transformed hex type lattice + CoverA_trans, & !< c/a ratio for transformed hex type lattice a_fcc, & !< lattice parameter a for fcc austenite a_bcc !< lattice paramater a for bcc martensite - mainProcess: if (worldrank == 0) then + mainProcess: if (worldrank == 0) then write(6,'(/,a)') ' <<<+- lattice init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" @@ -1373,7 +1375,7 @@ subroutine lattice_init if (iand(debug_level(debug_lattice),debug_levelBasic) /= 0_pInt) then write(6,'(a16,1x,i5)') ' # phases:',Nphases endif - + allocate(lattice_structure(Nphases),source = LATTICE_undefined_ID) allocate(trans_lattice_structure(Nphases),source = LATTICE_undefined_ID) allocate(lattice_C66(6,6,Nphases), source=0.0_pReal) @@ -1461,7 +1463,7 @@ subroutine lattice_init if (IO_isBlank(line)) cycle ! skip empty lines if (IO_getTag(line,'<','>') /= '') then ! stop at next part line = IO_read(fileUnit, .true.) ! reset IO_read - exit + exit endif if (IO_getTag(line,'[',']') /= '') then ! next section section = section + 1_pInt @@ -1615,7 +1617,7 @@ subroutine lattice_init end select endif enddo - + do i = 1_pInt,Nphases if ((CoverA(i) < 1.0_pReal .or. CoverA(i) > 2.0_pReal) & .and. lattice_structure(i) == LATTICE_hex_ID) call IO_error(131_pInt,el=i) ! checking physical significance of c/a @@ -1652,7 +1654,7 @@ subroutine lattice_initializeStructure(myPhase,CoverA,CoverA_trans,a_fcc,a_bcc) use IO, only: & IO_error, & IO_warning - + implicit none integer(pInt), intent(in) :: myPhase real(pReal), intent(in) :: & @@ -1689,7 +1691,7 @@ subroutine lattice_initializeStructure(myPhase,CoverA,CoverA_trans,a_fcc,a_bcc) lattice_C66(1:6,1:6,myPhase) = lattice_symmetrizeC66(lattice_structure(myPhase),& lattice_C66(1:6,1:6,myPhase)) - + lattice_mu(myPhase) = 0.2_pReal *( lattice_C66(1,1,myPhase) & - lattice_C66(1,2,myPhase) & + 3.0_pReal*lattice_C66(4,4,myPhase)) ! (C11iso-C12iso)/2 with C11iso=(3*C11+2*C12+4*C44)/5 and C12iso=(C11+4*C12-2*C44)/5 @@ -1735,7 +1737,7 @@ subroutine lattice_initializeStructure(myPhase,CoverA,CoverA_trans,a_fcc,a_bcc) lattice_trans_C66(1,3,myPhase) = c13bar lattice_trans_C66(3,3,myPhase) = c33bar lattice_trans_C66(4,4,myPhase) = c44bar - B - + lattice_trans_C66(1:6,1:6,myPhase) = lattice_symmetrizeC66(trans_lattice_structure(myPhase),& lattice_trans_C66(1:6,1:6,myPhase)) lattice_trans_mu(myPhase) = 0.2_pReal *( lattice_trans_C66(1,1,myPhase) & @@ -1772,7 +1774,7 @@ subroutine lattice_initializeStructure(myPhase,CoverA,CoverA_trans,a_fcc,a_bcc) lattice_hydrogenfluxDiffusion33(1:3,1:3,myPhase)) lattice_hydrogenfluxMobility33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& lattice_hydrogenfluxMobility33(1:3,1:3,myPhase)) - + select case(lattice_structure(myPhase)) !-------------------------------------------------------------------------------------------------- ! fcc @@ -1784,7 +1786,7 @@ subroutine lattice_initializeStructure(myPhase,CoverA,CoverA_trans,a_fcc,a_bcc) do i = 1_pInt,myNslip ! assign slip system vectors sd(1:3,i) = lattice_fcc_systemSlip(1:3,i) sn(1:3,i) = lattice_fcc_systemSlip(4:6,i) - enddo + enddo do i = 1_pInt,myNtwin ! assign twin system vectors and shears td(1:3,i) = lattice_fcc_systemTwin(1:3,i) tn(1:3,i) = lattice_fcc_systemTwin(4:6,i) @@ -1794,11 +1796,11 @@ subroutine lattice_initializeStructure(myPhase,CoverA,CoverA_trans,a_fcc,a_bcc) cd(1:3,i) = lattice_fcc_systemCleavage(1:3,i)/norm2(lattice_fcc_systemCleavage(1:3,i)) cn(1:3,i) = lattice_fcc_systemCleavage(4:6,i)/norm2(lattice_fcc_systemCleavage(4:6,i)) ct(1:3,i) = math_crossproduct(cd(1:3,i),cn(1:3,i)) - enddo + enddo ! Phase transformation select case(trans_lattice_structure(myPhase)) - case (LATTICE_bcc_ID) ! fcc to bcc transformation + case (LATTICE_bcc_ID) ! fcc to bcc transformation do i = 1_pInt,myNtrans Rtr(1:3,1:3,i) = math_axisAngleToR(lattice_fccTobcc_systemTrans(1:3,i), & ! Pitsch rotation lattice_fccTobcc_systemTrans(4,i)*INRAD) @@ -1894,7 +1896,7 @@ subroutine lattice_initializeStructure(myPhase,CoverA,CoverA_trans,a_fcc,a_bcc) cd(1:3,i) = lattice_bcc_systemCleavage(1:3,i)/norm2(lattice_bcc_systemCleavage(1:3,i)) cn(1:3,i) = lattice_bcc_systemCleavage(4:6,i)/norm2(lattice_bcc_systemCleavage(4:6,i)) ct(1:3,i) = math_crossproduct(cd(1:3,i),cn(1:3,i)) - enddo + enddo lattice_NslipSystem(1:lattice_maxNslipFamily,myPhase) = lattice_bcc_NslipSystem lattice_NtwinSystem(1:lattice_maxNtwinFamily,myPhase) = lattice_bcc_NtwinSystem lattice_NtransSystem(1:lattice_maxNtransFamily,myPhase) = lattice_bcc_NtransSystem @@ -1904,7 +1906,7 @@ subroutine lattice_initializeStructure(myPhase,CoverA,CoverA_trans,a_fcc,a_bcc) lattice_interactionSlipTwin(1:myNslip,1:myNtwin,myPhase) = lattice_bcc_interactionSlipTwin lattice_interactionTwinSlip(1:myNtwin,1:myNslip,myPhase) = lattice_bcc_interactionTwinSlip lattice_interactionTwinTwin(1:myNtwin,1:myNtwin,myPhase) = lattice_bcc_interactionTwinTwin - + !-------------------------------------------------------------------------------------------------- ! hex (including conversion from miller-bravais (a1=a2=a3=c) to miller (a, b, c) indices) case (LATTICE_hex_ID) @@ -1912,7 +1914,7 @@ subroutine lattice_initializeStructure(myPhase,CoverA,CoverA_trans,a_fcc,a_bcc) myNtwin = lattice_hex_Ntwin myNtrans = lattice_hex_Ntrans myNcleavage = lattice_hex_Ncleavage - do i = 1_pInt,myNslip ! assign slip system vectors + do i = 1_pInt,myNslip ! assign slip system vectors sd(1,i) = lattice_hex_systemSlip(1,i)*1.5_pReal ! direction [uvtw]->[3u/2 (u+2v)*sqrt(3)/2 w*(c/a)] sd(2,i) = (lattice_hex_systemSlip(1,i)+2.0_pReal*lattice_hex_systemSlip(2,i))*& 0.5_pReal*sqrt(3.0_pReal) @@ -1920,7 +1922,7 @@ subroutine lattice_initializeStructure(myPhase,CoverA,CoverA_trans,a_fcc,a_bcc) sn(1,i) = lattice_hex_systemSlip(5,i) ! plane (hkil)->(h (h+2k)/sqrt(3) l/(c/a)) sn(2,i) = (lattice_hex_systemSlip(5,i)+2.0_pReal*lattice_hex_systemSlip(6,i))/sqrt(3.0_pReal) sn(3,i) = lattice_hex_systemSlip(8,i)/CoverA - enddo + enddo do i = 1_pInt,myNtwin ! assign twin system vectors and shears td(1,i) = lattice_hex_systemTwin(1,i)*1.5_pReal td(2,i) = (lattice_hex_systemTwin(1,i)+2.0_pReal*lattice_hex_systemTwin(2,i))*& @@ -1940,7 +1942,7 @@ subroutine lattice_initializeStructure(myPhase,CoverA,CoverA_trans,a_fcc,a_bcc) ts(i) = 2.0_pReal*(CoverA*CoverA-2.0_pReal)/3.0_pReal/CoverA end select enddo - do i = 1_pInt, myNcleavage ! cleavage system vectors + do i = 1_pInt, myNcleavage ! cleavage system vectors cd(1,i) = lattice_hex_systemCleavage(1,i)*1.5_pReal ! direction [uvtw]->[3u/2 (u+2v)*sqrt(3)/2 w*(c/a)] cd(2,i) = (lattice_hex_systemCleavage(1,i)+2.0_pReal*lattice_hex_systemCleavage(2,i))*& 0.5_pReal*sqrt(3.0_pReal) @@ -1951,7 +1953,7 @@ subroutine lattice_initializeStructure(myPhase,CoverA,CoverA_trans,a_fcc,a_bcc) cn(3,i) = lattice_hex_systemCleavage(8,i)/CoverA cn(1:3,1) = cn(1:3,i)/norm2(cn(1:3,i)) ct(1:3,i) = math_crossproduct(cd(1:3,i),cn(1:3,i)) - enddo + enddo lattice_NslipSystem(1:lattice_maxNslipFamily,myPhase) = lattice_hex_NslipSystem lattice_NtwinSystem(1:lattice_maxNtwinFamily,myPhase) = lattice_hex_NtwinSystem lattice_NtransSystem(1:lattice_maxNtransFamily,myPhase) = lattice_hex_NtransSystem @@ -1994,7 +1996,7 @@ subroutine lattice_initializeStructure(myPhase,CoverA,CoverA_trans,a_fcc,a_bcc) cd(1:3,i) = lattice_iso_systemCleavage(1:3,i)/norm2(LATTICE_ortho_systemCleavage(1:3,i)) cn(1:3,i) = lattice_iso_systemCleavage(4:6,i)/norm2(LATTICE_ortho_systemCleavage(4:6,i)) ct(1:3,i) = math_crossproduct(cd(1:3,i),cn(1:3,i)) - enddo + enddo lattice_NcleavageSystem(1:lattice_maxNcleavageFamily,myPhase) = lattice_iso_NcleavageSystem !-------------------------------------------------------------------------------------------------- @@ -2008,7 +2010,7 @@ subroutine lattice_initializeStructure(myPhase,CoverA,CoverA_trans,a_fcc,a_bcc) cd(1:3,i) = lattice_iso_systemCleavage(1:3,i)/norm2(lattice_iso_systemCleavage(1:3,i)) cn(1:3,i) = lattice_iso_systemCleavage(4:6,i)/norm2(lattice_iso_systemCleavage(4:6,i)) ct(1:3,i) = math_crossproduct(cd(1:3,i),cn(1:3,i)) - enddo + enddo lattice_NcleavageSystem(1:lattice_maxNcleavageFamily,myPhase) = lattice_iso_NcleavageSystem !-------------------------------------------------------------------------------------------------- @@ -2028,7 +2030,7 @@ subroutine lattice_initializeStructure(myPhase,CoverA,CoverA_trans,a_fcc,a_bcc) do j = 1_pInt,lattice_NnonSchmid(myPhase) lattice_Sslip(1:3,1:3,2*j ,i,myPhase) = sns(1:3,1:3,1,j,i) lattice_Sslip(1:3,1:3,2*j+1,i,myPhase) = sns(1:3,1:3,2,j,i) - enddo + enddo do j = 1_pInt,1_pInt+2_pInt*lattice_NnonSchmid(myPhase) lattice_Sslip_v(1:6,j,i,myPhase) = & math_Mandel33to6(math_symmetric33(lattice_Sslip(1:3,1:3,j,i,myPhase))) @@ -2064,7 +2066,7 @@ subroutine lattice_initializeStructure(myPhase,CoverA,CoverA_trans,a_fcc,a_bcc) math_Mandel33to6(math_symmetric33(lattice_Scleavage(1:3,1:3,j,i,myPhase))) enddo enddo - + end subroutine lattice_initializeStructure @@ -2080,7 +2082,7 @@ pure function lattice_symmetrizeC66(struct,C66) integer(pInt) :: j,k lattice_symmetrizeC66 = 0.0_pReal - + select case(struct) case (LATTICE_iso_ID) forall(k=1_pInt:3_pInt) @@ -2093,7 +2095,7 @@ pure function lattice_symmetrizeC66(struct,C66) forall(j=1_pInt:3_pInt) lattice_symmetrizeC66(k,j) = C66(1,2) lattice_symmetrizeC66(k,k) = C66(1,1) lattice_symmetrizeC66(k+3_pInt,k+3_pInt) = C66(4,4) - end forall + end forall case (LATTICE_hex_ID) lattice_symmetrizeC66(1,1) = C66(1,1) lattice_symmetrizeC66(2,2) = C66(1,1) @@ -2136,7 +2138,7 @@ pure function lattice_symmetrizeC66(struct,C66) case default lattice_symmetrizeC66 = C66 end select - + end function lattice_symmetrizeC66 !-------------------------------------------------------------------------------------------------- @@ -2151,7 +2153,7 @@ pure function lattice_symmetrize33(struct,T33) integer(pInt) :: k lattice_symmetrize33 = 0.0_pReal - + select case(struct) case (LATTICE_iso_ID,LATTICE_fcc_ID,LATTICE_bcc_ID) forall(k=1_pInt:3_pInt) lattice_symmetrize33(k,k) = T33(1,1) @@ -2166,7 +2168,7 @@ pure function lattice_symmetrize33(struct,T33) case default lattice_symmetrize33 = T33 end select - + end function lattice_symmetrize33 @@ -2250,7 +2252,7 @@ pure function lattice_qDisorientation(Q1, Q2, struct) dQ = math_qMul(math_qConj(Q1),Q2) lattice_qDisorientation = dQ - select case(symmetry) + select case(symmetry) case (1_pInt,2_pInt) s = sum(lattice_NsymOperations(1:symmetry-1_pInt)) diff --git a/src/plastic_nonlocal.f90 b/src/plastic_nonlocal.f90 index b699c57ed..98e6af226 100644 --- a/src/plastic_nonlocal.f90 +++ b/src/plastic_nonlocal.f90 @@ -8,7 +8,7 @@ module plastic_nonlocal use prec, only: & pReal, & pInt - + implicit none private character(len=22), dimension(11), parameter, private :: & @@ -23,20 +23,20 @@ module plastic_nonlocal 'rhoDipEdge ', & 'rhoDipScrew ', & 'accumulatedshear ' ] !< list of "basic" microstructural state variables that are independent from other state variables - + character(len=16), dimension(3), parameter, private :: & - DEPENDENTSTATES = ['rhoForest ', & - 'tauThreshold ', & + DEPENDENTSTATES = ['rhoForest ', & + 'tauThreshold ', & 'tauBack ' ] !< list of microstructural state variables that depend on other state variables - + character(len=20), dimension(6), parameter, private :: & - OTHERSTATES = ['velocityEdgePos ', & - 'velocityEdgeNeg ', & + OTHERSTATES = ['velocityEdgePos ', & + 'velocityEdgeNeg ', & 'velocityScrewPos ', & 'velocityScrewNeg ', & 'maxDipoleHeightEdge ', & 'maxDipoleHeightScrew' ] !< list of other dependent state variables that are not updated by microstructure - + real(pReal), parameter, private :: & KB = 1.38e-23_pReal !< Physical parameter, Boltzmann constant in J/Kelvin @@ -48,13 +48,13 @@ module plastic_nonlocal integer(pInt), dimension(:,:), allocatable, target, public :: & plastic_nonlocal_sizePostResult !< size of each post result output - + character(len=64), dimension(:,:), allocatable, target, public :: & plastic_nonlocal_output !< name of each post result output - + integer(pInt), dimension(:), allocatable, target, public :: & - plastic_nonlocal_Noutput !< number of outputs per instance of this plasticity - + plastic_nonlocal_Noutput !< number of outputs per instance of this plasticity + integer(pInt), dimension(:,:), allocatable, private :: & iGamma, & !< state indices for accumulated shear iRhoF, & !< state indices for forest density @@ -66,16 +66,16 @@ module plastic_nonlocal iRhoD, & !< state indices for dipole density iV, & !< state indices for dislcation velocities iD !< state indices for stable dipole height - + integer(pInt), dimension(:), allocatable, public, protected :: & totalNslip !< total number of active slip systems for each instance - + integer(pInt), dimension(:,:), allocatable, private :: & Nslip, & !< number of active slip systems for each family and instance slipFamily, & !< lookup table relating active slip system to slip family for each instance slipSystemLattice, & !< lookup table relating active slip system index to lattice slip system index for each instance colinearSystem !< colinear system to the active slip system (only valid for fcc!) - + real(pReal), dimension(:), allocatable, private :: & atomicVolume, & !< atomic volume Dsd0, & !< prefactor for self-diffusion coefficient @@ -102,7 +102,7 @@ module plastic_nonlocal rhoSglRandomBinning, & linetensionEffect, & edgeJogFactor - + real(pReal), dimension(:,:), allocatable, private :: & rhoSglEdgePos0, & !< initial edge_pos dislocation density per slip system for each family and instance rhoSglEdgeNeg0, & !< initial edge_neg dislocation density per slip system for each family and instance @@ -115,40 +115,40 @@ module plastic_nonlocal burgersPerSlipFamily, & !< absolute length of burgers vector [m] for each family and instance burgers, & !< absolute length of burgers vector [m] for each slip system and instance interactionSlipSlip !< coefficients for slip-slip interaction for each interaction type and instance - + real(pReal), dimension(:,:,:), allocatable, private :: & minDipoleHeightPerSlipFamily, & !< minimum stable edge/screw dipole height for each family and instance minDipoleHeight, & !< minimum stable edge/screw dipole height for each slip system and instance - peierlsStressPerSlipFamily, & !< Peierls stress (edge and screw) - peierlsStress, & !< Peierls stress (edge and screw) + peierlsStressPerSlipFamily, & !< Peierls stress (edge and screw) + peierlsStress, & !< Peierls stress (edge and screw) forestProjectionEdge, & !< matrix of forest projections of edge dislocations for each instance forestProjectionScrew, & !< matrix of forest projections of screw dislocations for each instance interactionMatrixSlipSlip !< interaction matrix of the different slip systems for each instance - + real(pReal), dimension(:,:,:,:), allocatable, private :: & lattice2slip, & !< orthogonal transformation matrix from lattice coordinate system to slip coordinate system (passive rotation !!!) rhoDotEdgeJogsOutput, & sourceProbability - + real(pReal), dimension(:,:,:,:,:), allocatable, private :: & - rhoDotFluxOutput, & + rhoDotFluxOutput, & rhoDotMultiplicationOutput, & rhoDotSingle2DipoleGlideOutput, & rhoDotAthermalAnnihilationOutput, & rhoDotThermalAnnihilationOutput, & nonSchmidProjection !< combined projection of Schmid and non-Schmid contributions to the resolved shear stress (only for screws) - + real(pReal), dimension(:,:,:,:,:,:), allocatable, private :: & compatibility !< slip system compatibility between me and my neighbors - + real(pReal), dimension(:,:), allocatable, private :: & nonSchmidCoeff - + logical, dimension(:), allocatable, private :: & shortRangeStressCorrection, & !< flag indicating the use of the short range stress correction by a excess density gradient term probabilisticMultiplication - - enum, bind(c) + + enum, bind(c) enumerator :: undefined_ID, & rho_ID, & delta_ID, & @@ -235,9 +235,9 @@ module plastic_nonlocal accumulatedshear_ID, & dislocationstress_ID end enum - integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: & + integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: & plastic_nonlocal_outputID !< ID of each post result output - + public :: & plastic_nonlocal_init, & plastic_nonlocal_stateInit, & @@ -248,11 +248,11 @@ module plastic_nonlocal plastic_nonlocal_deltaState, & plastic_nonlocal_updateCompatibility, & plastic_nonlocal_postResults - + private :: & plastic_nonlocal_kinetics, & plastic_nonlocal_dislocationstress - + contains @@ -262,8 +262,8 @@ contains !-------------------------------------------------------------------------------------------------- subroutine plastic_nonlocal_init(fileUnit) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) -use math, only: math_Mandel3333to66, & - math_Voigt66to3333, & +use math, only: math_Mandel3333to66, & + math_Voigt66to3333, & math_mul3x3, & math_transpose33 use IO, only: IO_read, & @@ -330,9 +330,9 @@ integer(pInt) :: phase, & integer(pInt) :: sizeState, sizeDotState,sizeDependentState, sizeDeltaState - integer(pInt) :: NofMyPhase - - mainProcess: if (worldrank == 0) then + integer(pInt) :: NofMyPhase + + mainProcess: if (worldrank == 0) then write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_NONLOCAL_label//' init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" @@ -406,11 +406,11 @@ allocate(nonSchmidCoeff(lattice_maxNnonSchmid,maxNinstances), s do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to line = IO_read(fileUnit) enddo - + parsingFile: do while (trim(line) /= IO_EOF) ! read through phases of phase part line = IO_read(fileUnit) if (IO_isBlank(line)) cycle ! skip empty lines - if (IO_getTag(line,'<','>') /= '') then ! stop at next part + if (IO_getTag(line,'<','>') /= '') then ! stop at next part line = IO_read(fileUnit, .true.) ! reset IO_read exit endif @@ -978,7 +978,7 @@ allocate(nonSchmidCoeff(lattice_maxNnonSchmid,maxNinstances), s sanityChecks: do phase = 1_pInt, size(phase_plasticity) myPhase: if (phase_plasticity(phase) == PLASTICITY_NONLOCAL_ID) then - instance = phase_plasticityInstance(phase) + instance = phase_plasticityInstance(phase) if (sum(Nslip(:,instance)) <= 0_pInt) & call IO_error(211_pInt,ext_msg='Nslip ('//PLASTICITY_NONLOCAL_label//')') do o = 1_pInt,maxval(phase_Noutput) @@ -1065,13 +1065,13 @@ allocate(nonSchmidCoeff(lattice_maxNnonSchmid,maxNinstances), s call IO_error(211_pInt,ext_msg='CFLfactor ('//PLASTICITY_NONLOCAL_label//')') if (fEdgeMultiplication(instance) < 0.0_pReal .or. fEdgeMultiplication(instance) > 1.0_pReal) & call IO_error(211_pInt,ext_msg='edgemultiplicationfactor ('//PLASTICITY_NONLOCAL_label//')') - - + + !*** determine total number of active slip systems Nslip(1:lattice_maxNslipFamily,instance) = min(lattice_NslipSystem(1:lattice_maxNslipFamily,phase), & Nslip(1:lattice_maxNslipFamily,instance) ) ! we can't use more slip systems per family than specified in lattice totalNslip(instance) = sum(Nslip(1:lattice_maxNslipFamily,instance)) - endif myPhase + endif myPhase enddo sanityChecks @@ -1116,13 +1116,13 @@ allocate(compatibility(2,maxTotalNslip,maxTotalNslip,mesh_maxNipNeighbors,mesh_m allocate(peierlsStress(maxTotalNslip,2,maxNinstances), source=0.0_pReal) allocate(colinearSystem(maxTotalNslip,maxNinstances), source=0_pInt) allocate(nonSchmidProjection(3,3,4,maxTotalNslip,maxNinstances), source=0.0_pReal) - + initializeInstances: do phase = 1_pInt, size(phase_plasticity) NofMyPhase=count(material_phase==phase) myPhase2: if (phase_plasticity(phase) == PLASTICITY_NONLOCAL_ID .and. NofMyPhase/=0) then instance = phase_plasticityInstance(phase) !*** Inverse lookup of my slip system family and the slip system in lattice - + l = 0_pInt do f = 1_pInt,lattice_maxNslipFamily do s = 1_pInt,Nslip(f,instance) @@ -1130,10 +1130,10 @@ allocate(nonSchmidProjection(3,3,4,maxTotalNslip,maxNinstances), slipFamily(l,instance) = f slipSystemLattice(l,instance) = sum(lattice_NslipSystem(1:f-1_pInt, phase)) + s enddo; enddo - - + + !*** determine size of state array - + ns = totalNslip(instance) sizeDotState = int(size(BASICSTATES),pInt) * ns @@ -1193,10 +1193,10 @@ allocate(nonSchmidProjection(3,3,4,maxTotalNslip,maxNinstances), enddo if (iD(ns,2,instance) /= sizeState) & ! check if last index is equal to size of state call IO_error(0_pInt, ext_msg = 'state indices not properly set ('//PLASTICITY_NONLOCAL_label//')') - - + + !*** determine size of postResults array - + outputsLoop: do o = 1_pInt,plastic_nonlocal_Noutput(instance) select case(plastic_nonlocal_outputID(o,instance)) case( rho_ID, & @@ -1287,13 +1287,13 @@ allocate(nonSchmidProjection(3,3,4,maxTotalNslip,maxNinstances), mySize = 6_pInt case default end select - - if (mySize > 0_pInt) then ! any meaningful output found + + if (mySize > 0_pInt) then ! any meaningful output found plastic_nonlocal_sizePostResult(o,instance) = mySize plastic_nonlocal_sizePostResults(instance) = plastic_nonlocal_sizePostResults(instance) + mySize endif enddo outputsLoop - + plasticState(phase)%sizeState = sizeState plasticState(phase)%sizeDotState = sizeDotState plasticState(phase)%sizeDeltaState = sizeDeltaState @@ -1326,63 +1326,63 @@ allocate(nonSchmidProjection(3,3,4,maxTotalNslip,maxNinstances), plasticState(phase)%dotState(iGamma(1,instance):iGamma(ns,instance),1:NofMyPhase) plasticState(phase)%accumulatedSlip => & plasticState(phase)%state (iGamma(1,instance):iGamma(ns,instance),1:NofMyPhase) - - do s1 = 1_pInt,ns + + do s1 = 1_pInt,ns f = slipFamily(s1,instance) - + !*** burgers vector, mean free path prefactor and minimum dipole distance for each slip system - + burgers(s1,instance) = burgersPerSlipFamily(f,instance) lambda0(s1,instance) = lambda0PerSlipFamily(f,instance) minDipoleHeight(s1,1:2,instance) = minDipoleHeightPerSlipFamily(f,1:2,instance) peierlsStress(s1,1:2,instance) = peierlsStressPerSlipFamily(f,1:2,instance) - + do s2 = 1_pInt,ns - + !*** calculation of forest projections for edge and screw dislocations. s2 acts as forest for s1 - + forestProjectionEdge(s1,s2,instance) & = abs(math_mul3x3(lattice_sn(1:3,slipSystemLattice(s1,instance),phase), & lattice_st(1:3,slipSystemLattice(s2,instance),phase))) ! forest projection of edge dislocations is the projection of (t = b x n) onto the slip normal of the respective slip plane - + forestProjectionScrew(s1,s2,instance) & = abs(math_mul3x3(lattice_sn(1:3,slipSystemLattice(s1,instance),phase), & lattice_sd(1:3,slipSystemLattice(s2,instance),phase))) ! forest projection of screw dislocations is the projection of b onto the slip normal of the respective splip plane - + !*** calculation of interaction matrices - + interactionMatrixSlipSlip(s1,s2,instance) & = interactionSlipSlip(lattice_interactionSlipSlip(slipSystemLattice(s1,instance), & slipSystemLattice(s2,instance), & phase), instance) - + !*** colinear slip system (only makes sense for fcc like it is defined here) - + if (lattice_interactionSlipSlip(slipSystemLattice(s1,instance), & slipSystemLattice(s2,instance), & phase) == 3_pInt) then colinearSystem(s1,instance) = s2 endif - + enddo - + !*** rotation matrix from lattice configuration to slip system - + lattice2slip(1:3,1:3,s1,instance) & = math_transpose33( reshape([ lattice_sd(1:3, slipSystemLattice(s1,instance), phase), & -lattice_st(1:3, slipSystemLattice(s1,instance), phase), & lattice_sn(1:3, slipSystemLattice(s1,instance), phase)], [3,3])) enddo - - + + !*** combined projection of Schmid and non-Schmid contributions to the resolved shear stress (only for screws) - !* four types t: + !* four types t: !* 1) positive screw at positive resolved stress !* 2) positive screw at negative resolved stress !* 3) negative screw at positive resolved stress !* 4) negative screw at negative resolved stress - - do s = 1_pInt,ns + + do s = 1_pInt,ns do l = 1_pInt,lattice_NnonSchmid(phase) nonSchmidProjection(1:3,1:3,1,s,instance) = nonSchmidProjection(1:3,1:3,1,s,instance) & + nonSchmidCoeff(l,instance) * lattice_Sslip(1:3,1:3,2*l,slipSystemLattice(s,instance),phase) @@ -1398,7 +1398,7 @@ allocate(nonSchmidProjection(3,3,4,maxTotalNslip,maxNinstances), call plastic_nonlocal_aTolState(phase,instance) endif myPhase2 - + enddo initializeInstances end subroutine plastic_nonlocal_init @@ -1426,7 +1426,7 @@ implicit none integer(pInt) :: e, & i, & - ns, & ! short notation for total number of active slip systems + ns, & ! short notation for total number of active slip systems f, & ! index of lattice family from, & upto, & @@ -1436,7 +1436,7 @@ integer(pInt) :: e, & instance, & maxNinstances real(pReal), dimension(2) :: noise -real(pReal), dimension(4) :: rnd +real(pReal), dimension(4) :: rnd real(pReal) meanDensity, & totalVolume, & densityBinning, & @@ -1447,7 +1447,7 @@ maxNinstances = int(count(phase_plasticity == PLASTICITY_NONLOCAL_ID),pInt) do instance = 1_pInt,maxNinstances ns = totalNslip(instance) - ! randomly distribute dislocation segments on random slip system and of random type in the volume + ! randomly distribute dislocation segments on random slip system and of random type in the volume if (rhoSglRandom(instance) > 0.0_pReal) then ! get the total volume of the instance @@ -1526,7 +1526,7 @@ subroutine plastic_nonlocal_aTolState(ph,instance) plasticState implicit none - integer(pInt), intent(in) :: & + integer(pInt), intent(in) :: & instance, & !< number specifying the instance of the plasticity ph integer(pInt) :: & @@ -1540,7 +1540,7 @@ subroutine plastic_nonlocal_aTolState(ph,instance) end forall forall (c = 1_pInt:2_pInt) & plasticState(ph)%aTolState(iRhoD(1:ns,c,instance)) = aTolRho(instance) - + plasticState(ph)%aTolState(iGamma(1:ns,instance)) = aTolShear(instance) end subroutine plastic_nonlocal_aTolState @@ -1606,8 +1606,8 @@ real(pReal), dimension(3,3), intent(in) :: & integer(pInt) neighbor_el, & ! element number of neighboring material point neighbor_ip, & ! integration point of neighboring material point - instance, & ! my instance of this plasticity - neighbor_instance, & ! instance of this plasticity of neighboring material point + instance, & ! my instance of this plasticity + neighbor_instance, & ! instance of this plasticity of neighboring material point neighbor_phase, & ns, & ! total number of active slip systems at my material point neighbor_ns, & ! total number of active slip systems at neighboring material point @@ -1678,25 +1678,25 @@ where (abs(rhoDip) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(instance forall (s = 1_pInt:ns) & rhoForest(s) = dot_product((sum(abs(rhoSgl(1:ns,[1,2,5,6])),2) + rhoDip(1:ns,1)), & - forestProjectionEdge(s,1:ns,instance)) & + forestProjectionEdge(s,1:ns,instance)) & + dot_product((sum(abs(rhoSgl(1:ns,[3,4,7,8])),2) + rhoDip(1:ns,2)), & forestProjectionScrew(s,1:ns,instance)) -!*** calculate the threshold shear stress for dislocation slip -!*** coefficients are corrected for the line tension effect +!*** calculate the threshold shear stress for dislocation slip +!*** coefficients are corrected for the line tension effect !*** (see Kubin,Devincre,Hoc; 2008; Modeling dislocation storage rates and mean free paths in face-centered cubic crystals) myInteractionMatrix = 0.0_pReal myInteractionMatrix(1:ns,1:ns) = interactionMatrixSlipSlip(1:ns,1:ns,instance) if (lattice_structure(ph) == LATTICE_bcc_ID .or. lattice_structure(ph) == LATTICE_fcc_ID) then ! only fcc and bcc - do s = 1_pInt,ns + do s = 1_pInt,ns myRhoForest = max(rhoForest(s),significantRho(instance)) correction = ( 1.0_pReal - linetensionEffect(instance) & + linetensionEffect(instance) & * log(0.35_pReal * burgers(s,instance) * sqrt(myRhoForest)) & / log(0.35_pReal * burgers(s,instance) * 1e6_pReal)) ** 2.0_pReal - myInteractionMatrix(s,1:ns) = correction * myInteractionMatrix(s,1:ns) + myInteractionMatrix(s,1:ns) = correction * myInteractionMatrix(s,1:ns) enddo endif forall (s = 1_pInt:ns) & @@ -1715,9 +1715,9 @@ if (.not. phase_localPlasticity(ph) .and. shortRangeStressCorrection(instance)) rhoExcess(1,1:ns) = rhoSgl(1:ns,1) - rhoSgl(1:ns,2) rhoExcess(2,1:ns) = rhoSgl(1:ns,3) - rhoSgl(1:ns,4) FVsize = mesh_ipVolume(ip,el) ** (1.0_pReal/3.0_pReal) - + !* loop through my neighborhood and get the connection vectors (in lattice frame) and the excess densities - + nRealNeighbors = 0_pInt neighbor_rhoTotal = 0.0_pReal do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el)))) @@ -1769,17 +1769,17 @@ if (.not. phase_localPlasticity(ph) .and. shortRangeStressCorrection(instance)) neighbor_rhoExcess(1:2,1:ns,n) = rhoExcess endif enddo - + !* loop through the slip systems and calculate the dislocation gradient by !* 1. interpolation of the excess density in the neighorhood !* 2. interpolation of the dead dislocation density in the central volume - + m(1:3,1:ns,1) = lattice_sd(1:3,slipSystemLattice(1:ns,instance),ph) m(1:3,1:ns,2) = -lattice_st(1:3,slipSystemLattice(1:ns,instance),ph) do s = 1_pInt,ns - + !* gradient from interpolation of neighboring excess density do c = 1_pInt,2_pInt @@ -1797,23 +1797,23 @@ if (.not. phase_localPlasticity(ph) .and. shortRangeStressCorrection(instance)) rhoExcessGradient(c) = math_mul3x3(m(1:3,s,c), & math_mul33x3(invConnections,rhoExcessDifferences)) enddo - + !* plus gradient from deads - + do t = 1_pInt,4_pInt c = (t - 1_pInt) / 2_pInt + 1_pInt rhoExcessGradient(c) = rhoExcessGradient(c) + rhoSgl(s,t+4_pInt) / FVsize enddo !* normalized with the total density - + rhoExcessGradient_over_rho = 0.0_pReal forall (c = 1_pInt:2_pInt) & rhoTotal(c) = (sum(abs(rhoSgl(s,[2*c-1,2*c,2*c+3,2*c+4]))) + rhoDip(s,c) & + sum(neighbor_rhoTotal(c,s,:))) / real(1_pInt + nRealNeighbors,pReal) forall (c = 1_pInt:2_pInt, rhoTotal(c) > 0.0_pReal) & rhoExcessGradient_over_rho(c) = rhoExcessGradient(c) / rhoTotal(c) - + !* gives the local stress correction when multiplied with a factor tauBack(s) = - lattice_mu(ph) * burgers(s,instance) / (2.0_pReal * pi) & @@ -1844,7 +1844,7 @@ end subroutine plastic_nonlocal_microstructure !-------------------------------------------------------------------------------------------------- -!> @brief calculates kinetics +!> @brief calculates kinetics !-------------------------------------------------------------------------------------------------- subroutine plastic_nonlocal_kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, & tauThreshold, c, Temperature, ip, el) @@ -1880,7 +1880,7 @@ real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1_pInt integer(pInt) :: instance, & !< current instance of this plasticity ns, & !< short notation for the total number of active slip systems s !< index of my current slip system -real(pReal) tauRel_P, & +real(pReal) tauRel_P, & tauRel_S, & tauEff, & !< effective shear stress tPeierls, & !< waiting time in front of a peierls barriers @@ -1918,7 +1918,7 @@ if (Temperature > 0.0_pReal) then !* Peierls contribution !* Effective stress includes non Schmid constributions !* The derivative only gives absolute values; the correct sign is taken care of in the formula for the derivative of the velocity - + tauEff = max(0.0_pReal, abs(tauNS(s)) - tauThreshold(s)) ! ensure that the effective stress is positive meanfreepath_P = burgers(s,instance) jumpWidth_P = burgers(s,instance) @@ -1933,7 +1933,7 @@ if (Temperature > 0.0_pReal) then if (tauEff < criticalStress_P) then dtPeierls_dtau = tPeierls * pParam(instance) * qParam(instance) * activationVolume_P / (KB * Temperature) & * (1.0_pReal - tauRel_P**pParam(instance))**(qParam(instance)-1.0_pReal) & - * tauRel_P**(pParam(instance)-1.0_pReal) + * tauRel_P**(pParam(instance)-1.0_pReal) else dtPeierls_dtau = 0.0_pReal endif @@ -1957,32 +1957,32 @@ if (Temperature > 0.0_pReal) then dtSolidSolution_dtau = tSolidSolution * pParam(instance) * qParam(instance) & * activationVolume_S / (KB * Temperature) & * (1.0_pReal - tauRel_S**pParam(instance))**(qParam(instance)-1.0_pReal) & - * tauRel_S**(pParam(instance)-1.0_pReal) + * tauRel_S**(pParam(instance)-1.0_pReal) else dtSolidSolution_dtau = 0.0_pReal endif !* viscous glide velocity - + tauEff = abs(tau(s)) - tauThreshold(s) mobility = burgers(s,instance) / viscosity(instance) vViscous = mobility * tauEff - !* Mean velocity results from waiting time at peierls barriers and solid solution obstacles with respective meanfreepath of - !* free flight at glide velocity in between. + !* Mean velocity results from waiting time at peierls barriers and solid solution obstacles with respective meanfreepath of + !* free flight at glide velocity in between. !* adopt sign from resolved stress v(s) = sign(1.0_pReal,tau(s)) & / (tPeierls / meanfreepath_P + tSolidSolution / meanfreepath_S + 1.0_pReal / vViscous) dv_dtau(s) = v(s) * v(s) * (dtSolidSolution_dtau / meanfreepath_S & + mobility / (vViscous * vViscous)) - dv_dtauNS(s) = v(s) * v(s) * dtPeierls_dtau / meanfreepath_P + dv_dtauNS(s) = v(s) * v(s) * dtPeierls_dtau / meanfreepath_P endif enddo endif - + #ifndef _OPENMP if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt & @@ -2051,7 +2051,7 @@ integer(pInt) instance, & sLattice !< index of my current slip system according to lattice order real(pReal), dimension(3,3,3,3) :: dLp_dTstar3333 !< derivative of Lp with respect to Tstar (3x3x3x3 matrix) real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1_pInt,ip,el))),8) :: & - rhoSgl !< single dislocation densities (including blocked) + rhoSgl !< single dislocation densities (including blocked) real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1_pInt,ip,el))),4) :: & v, & !< velocity tauNS, & !< resolved shear stress including non Schmid and backstress terms @@ -2075,7 +2075,7 @@ instance = phase_plasticityInstance(ph) ns = totalNslip(instance) -!*** shortcut to state variables +!*** shortcut to state variables forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) @@ -2114,7 +2114,7 @@ tau = tau + tauBack !*** get dislocation velocity and its tangent and store the velocity in the state array -! edges +! edges call plastic_nonlocal_kinetics(v(1:ns,1), dv_dtau(1:ns,1), dv_dtauNS(1:ns,1), & tau(1:ns), tauNS(1:ns,1), tauThreshold(1:ns), & 1_pInt, Temperature, ip, el) @@ -2153,7 +2153,7 @@ forall (s = 1_pInt:ns, t = 5_pInt:8_pInt, rhoSgl(s,t) * v(s,t-4_pInt) < 0.0_pRea gdotTotal = sum(rhoSgl(1:ns,1:4) * v, 2) * burgers(1:ns,instance) do s = 1_pInt,ns - sLattice = slipSystemLattice(s,instance) + sLattice = slipSystemLattice(s,instance) Lp = Lp + gdotTotal(s) * lattice_Sslip(1:3,1:3,1,sLattice,ph) ! Schmid contributions to tangent @@ -2167,14 +2167,14 @@ do s = 1_pInt,ns forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt,k=1_pInt:3_pInt,l=1_pInt:3_pInt) & dLp_dTstar3333(i,j,k,l) = dLp_dTstar3333(i,j,k,l) & + lattice_Sslip(i,j,1,sLattice,ph) & - * ( nonSchmidProjection(k,l,1,s,instance) * rhoSgl(s,3) * dv_dtauNS(s,3) & + * ( nonSchmidProjection(k,l,1,s,instance) * rhoSgl(s,3) * dv_dtauNS(s,3) & + nonSchmidProjection(k,l,3,s,instance) * rhoSgl(s,4) * dv_dtauNS(s,4) ) & * burgers(s,instance) else forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt,k=1_pInt:3_pInt,l=1_pInt:3_pInt) & dLp_dTstar3333(i,j,k,l) = dLp_dTstar3333(i,j,k,l) & + lattice_Sslip(i,j,1,sLattice,ph) & - * ( nonSchmidProjection(k,l,2,s,instance) * rhoSgl(s,3) * dv_dtauNS(s,3) & + * ( nonSchmidProjection(k,l,2,s,instance) * rhoSgl(s,3) * dv_dtauNS(s,3) & + nonSchmidProjection(k,l,4,s,instance) * rhoSgl(s,4) * dv_dtauNS(s,4) ) & * burgers(s,instance) endif @@ -2265,7 +2265,7 @@ real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,e ns = totalNslip(instance) -!*** shortcut to state variables +!*** shortcut to state variables forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) rhoSgl(s,t) = max(plasticState(ph)%state(iRhoU(s,t,instance),of), 0.0_pReal) ! ensure positive single mobile densities @@ -2311,7 +2311,7 @@ enddo !*** calculate limits for stable dipole height do s = 1_pInt,ns - sLattice = slipSystemLattice(s,instance) + sLattice = slipSystemLattice(s,instance) tau(s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,ph)) + tauBack(s) if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal enddo @@ -2324,7 +2324,7 @@ dUpper(1:ns,2) = lattice_mu(ph) * burgers(1:ns,instance) / (4.0_pReal * pi * abs forall (c = 1_pInt:2_pInt) where(sqrt(rhoSgl(1:ns,2*c-1)+rhoSgl(1:ns,2*c)+& abs(rhoSgl(1:ns,2*c+3))+abs(rhoSgl(1:ns,2*c+4))+rhoDip(1:ns,c)) >= tiny(0.0_pReal)) & - dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) & + dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) & + abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), & dUpper(1:ns,c)) end forall @@ -2347,7 +2347,7 @@ forall (t=1_pInt:4_pInt) & !*** store new maximum dipole height in state forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) & - plasticState(ph)%state(iD(s,c,instance),of) = dUpper(s,c) + plasticState(ph)%state(iD(s,c,instance),of) = dUpper(s,c) @@ -2400,7 +2400,7 @@ use math, only: math_mul6x6, & math_mul33x33, & math_inv33, & math_det33, & - math_transpose33, & + math_transpose33, & pi use mesh, only: mesh_NcpElems, & mesh_maxNips, & @@ -2426,7 +2426,7 @@ use lattice, only: lattice_Sslip_v, & lattice_mu, & lattice_nu, & lattice_structure, & - LATTICE_bcc_ID, & + LATTICE_bcc_ID, & LATTICE_fcc_ID implicit none @@ -2438,14 +2438,14 @@ real(pReal), intent(in) :: Temperature, & timestep !< substepped crystallite time increment real(pReal), dimension(6), intent(in) :: Tstar_v !< current 2nd Piola-Kirchhoff stress in Mandel notation real(pReal), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & - subfrac !< fraction of timestep at the beginning of the substepped crystallite time increment + subfrac !< fraction of timestep at the beginning of the substepped crystallite time increment real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & Fe, & !< elastic deformation gradient Fp !< plastic deformation gradient - + !*** local variables -integer(pInt) :: ph, & +integer(pInt) :: ph, & instance, & !< current instance of this plasticity neighbor_instance, & !< instance of my neighbor's plasticity ns, & !< short notation for the total number of active slip systems @@ -2494,7 +2494,7 @@ real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1_pInt nSources real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1_pInt,ip,el))),2) :: & rhoDip, & !< current dipole dislocation densities (screw and edge dipoles) - rhoDipOriginal, & + rhoDipOriginal, & dLower, & !< minimum stable dipole distance for edges and screws dUpper !< current maximum stable dipole distance for edges and screws real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(1_pInt,ip,el))),4) :: & @@ -2512,11 +2512,11 @@ real(pReal) area, & transmissivity, & !< overall transmissivity of dislocation flux to neighboring material point lineLength, & !< dislocation line length leaving the current interface selfDiffusion, & !< self diffusion - rnd, & + rnd, & meshlength logical considerEnteringFlux, & considerLeavingFlux - + p = phaseAt(1,ip,el) o = phasememberAt(1,ip,el) @@ -2538,7 +2538,7 @@ tau = 0.0_pReal gdot = 0.0_pReal -!*** shortcut to state variables +!*** shortcut to state variables forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) @@ -2572,7 +2572,7 @@ if (numerics_timeSyncing) then .or. abs(rhoSgl0) < significantRho(instance)) & rhoSgl0 = 0.0_pReal endif - + !*** sanity check for timestep @@ -2605,7 +2605,7 @@ forall (t = 1_pInt:4_pInt) & !*** calculate limits for stable dipole height do s = 1_pInt,ns ! loop over slip systems - sLattice = slipSystemLattice(s,instance) + sLattice = slipSystemLattice(s,instance) tau(s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,ph)) + tauBack(s) if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal enddo @@ -2618,7 +2618,7 @@ dUpper(1:ns,2) = lattice_mu(ph) * burgers(1:ns,instance) & forall (c = 1_pInt:2_pInt) where(sqrt(rhoSgl(1:ns,2*c-1)+rhoSgl(1:ns,2*c)+& abs(rhoSgl(1:ns,2*c+3))+abs(rhoSgl(1:ns,2*c+4))+rhoDip(1:ns,c)) >= tiny(0.0_pReal)) & - dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) & + dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) & + abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), & dUpper(1:ns,c)) end forall @@ -2692,7 +2692,7 @@ if (.not. phase_localPlasticity(material_phase(1_pInt,ip,el))) then .and. CFLfactor(instance) * abs(v) * timestep & > mesh_ipVolume(ip,el) / maxval(mesh_ipArea(:,ip,el)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here) #ifndef _OPENMP - if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt) then + if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt) then write(6,'(a,i5,a,i2)') '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip write(6,'(a,e10.3,a,e10.3)') '<< CONST >> velocity is at ', & maxval(abs(v), abs(gdot) > 0.0_pReal & @@ -2709,15 +2709,15 @@ if (.not. phase_localPlasticity(material_phase(1_pInt,ip,el))) then !*** be aware of the definition of lattice_st = lattice_sd x lattice_sn !!! !*** opposite sign to our p vector in the (s,p,n) triplet !!! - + m(1:3,1:ns,1) = lattice_sd(1:3, slipSystemLattice(1:ns,instance), ph) m(1:3,1:ns,2) = -lattice_sd(1:3, slipSystemLattice(1:ns,instance), ph) m(1:3,1:ns,3) = -lattice_st(1:3, slipSystemLattice(1:ns,instance), ph) m(1:3,1:ns,4) = lattice_st(1:3, slipSystemLattice(1:ns,instance), ph) - + my_Fe = Fe(1:3,1:3,1_pInt,ip,el) my_F = math_mul33x33(my_Fe, Fp(1:3,1:3,1_pInt,ip,el)) - + do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el)))) ! loop through my neighbors ! write(6,*) 'c' neighbor_el = mesh_ipNeighborhood(1,n,ip,el) @@ -2730,7 +2730,7 @@ if (.not. phase_localPlasticity(material_phase(1_pInt,ip,el))) then opposite_el = mesh_ipNeighborhood(1,opposite_neighbor,ip,el) opposite_ip = mesh_ipNeighborhood(2,opposite_neighbor,ip,el) opposite_n = mesh_ipNeighborhood(3,opposite_neighbor,ip,el) - + if (neighbor_n > 0_pInt) then ! if neighbor exists, average deformation gradient neighbor_instance = phase_plasticityInstance(material_phase(1_pInt,neighbor_ip,neighbor_el)) neighbor_Fe = Fe(1:3,1:3,1_pInt,neighbor_ip,neighbor_el) @@ -2739,26 +2739,26 @@ if (.not. phase_localPlasticity(material_phase(1_pInt,ip,el))) then else ! if no neighbor, take my value as average Favg = my_F endif - + !* FLUX FROM MY NEIGHBOR TO ME - !* This is only considered, if I have a neighbor of nonlocal plasticity - !* (also nonlocal constitutive law with local properties) that is at least a little bit + !* This is only considered, if I have a neighbor of nonlocal plasticity + !* (also nonlocal constitutive law with local properties) that is at least a little bit !* compatible. !* If it's not at all compatible, no flux is arriving, because everything is dammed in front of !* my neighbor's interface. - !* The entering flux from my neighbor will be distributed on my slip systems according to the + !* The entering flux from my neighbor will be distributed on my slip systems according to the !*compatibility - + considerEnteringFlux = .false. - neighbor_v = 0.0_pReal ! needed for check of sign change in flux density below + neighbor_v = 0.0_pReal ! needed for check of sign change in flux density below neighbor_rhoSgl = 0.0_pReal if (neighbor_n > 0_pInt) then if (phase_plasticity(material_phase(1,neighbor_ip,neighbor_el)) == PLASTICITY_NONLOCAL_ID & .and. any(compatibility(:,:,:,n,ip,el) > 0.0_pReal)) & considerEnteringFlux = .true. endif - + if (considerEnteringFlux) then if(numerics_timeSyncing .and. (subfrac(1_pInt,neighbor_ip,neighbor_el) /= subfrac(1_pInt,ip,el))) & then ! for timesyncing: in case of a timestep at the interface we have to use "state0" to make sure that fluxes n both sides are equal @@ -2790,7 +2790,7 @@ if (.not. phase_localPlasticity(material_phase(1_pInt,ip,el))) then c = (t + 1_pInt) / 2 topp = t + mod(t,2_pInt) - mod(t+1_pInt,2_pInt) if (neighbor_v(s,t) * math_mul3x3(m(1:3,s,t), normal_neighbor2me) > 0.0_pReal & ! flux from my neighbor to me == entering flux for me - .and. v(s,t) * neighbor_v(s,t) >= 0.0_pReal ) then ! ... only if no sign change in flux density + .and. v(s,t) * neighbor_v(s,t) >= 0.0_pReal ) then ! ... only if no sign change in flux density lineLength = neighbor_rhoSgl(s,t) * neighbor_v(s,t) & * math_mul3x3(m(1:3,s,t), normal_neighbor2me) * area ! positive line length that wants to enter through this interface where (compatibility(c,1_pInt:ns,s,n,ip,el) > 0.0_pReal) & ! positive compatibility... @@ -2805,16 +2805,16 @@ if (.not. phase_localPlasticity(material_phase(1_pInt,ip,el))) then enddo enddo endif - - + + !* FLUX FROM ME TO MY NEIGHBOR - !* This is not considered, if my opposite neighbor has a different constitutive law than nonlocal (still considered for nonlocal law with lcal properties). + !* This is not considered, if my opposite neighbor has a different constitutive law than nonlocal (still considered for nonlocal law with lcal properties). !* Then, we assume, that the opposite(!) neighbor sends an equal amount of dislocations to me. !* So the net flux in the direction of my neighbor is equal to zero: !* leaving flux to neighbor == entering flux from opposite neighbor !* In case of reduced transmissivity, part of the leaving flux is stored as dead dislocation density. !* That means for an interface of zero transmissivity the leaving flux is fully converted to dead dislocations. - + considerLeavingFlux = .true. if (opposite_n > 0_pInt) then if (phase_plasticity(material_phase(1,opposite_ip,opposite_el)) /= PLASTICITY_NONLOCAL_ID) & @@ -2822,10 +2822,10 @@ if (.not. phase_localPlasticity(material_phase(1_pInt,ip,el))) then endif if (considerLeavingFlux) then - - !* timeSyncing mode: If the central ip has zero subfraction, always use "state0". This is needed in case of + + !* timeSyncing mode: If the central ip has zero subfraction, always use "state0". This is needed in case of !* a synchronization step for the central ip, because then "state" contains the values at the end of the - !* previously converged full time step. Also, if either me or my neighbor has zero subfraction, we have to + !* previously converged full time step. Also, if either me or my neighbor has zero subfraction, we have to !* use "state0" to make sure that fluxes on both sides of the (potential) timestep are equal. my_rhoSgl = rhoSgl my_v = v @@ -2840,14 +2840,14 @@ if (.not. phase_localPlasticity(material_phase(1_pInt,ip,el))) then endif endif endif - + normal_me2neighbor_defConf = math_det33(Favg) & - * math_mul33x3(math_inv33(math_transpose33(Favg)), & + * math_mul33x3(math_inv33(math_transpose33(Favg)), & mesh_ipAreaNormal(1:3,n,ip,el)) ! calculate the normal of the interface in (average) deformed configuration (pointing from me to my neighbor!!!) normal_me2neighbor = math_mul33x3(math_transpose33(my_Fe), normal_me2neighbor_defConf) & / math_det33(my_Fe) ! interface normal in my lattice configuration area = mesh_ipArea(n,ip,el) * norm2(normal_me2neighbor) - normal_me2neighbor = normal_me2neighbor / norm2(normal_me2neighbor) ! normalize the surface normal to unit length + normal_me2neighbor = normal_me2neighbor / norm2(normal_me2neighbor) ! normalize the surface normal to unit length do s = 1_pInt,ns do t = 1_pInt,4_pInt c = (t + 1_pInt) / 2_pInt @@ -2866,9 +2866,9 @@ if (.not. phase_localPlasticity(material_phase(1_pInt,ip,el))) then endif enddo enddo - endif - - enddo ! neighbor loop + endif + + enddo ! neighbor loop endif @@ -2906,7 +2906,7 @@ enddo rhoDotAthermalAnnihilation = 0.0_pReal -forall (c=1_pInt:2_pInt) & +forall (c=1_pInt:2_pInt) & rhoDotAthermalAnnihilation(1:ns,c+8_pInt) = -2.0_pReal * dLower(1:ns,c) / burgers(1:ns,instance) & * ( 2.0_pReal * (rhoSgl(1:ns,2*c-1) * abs(gdot(1:ns,2*c)) + rhoSgl(1:ns,2*c) * abs(gdot(1:ns,2*c-1))) & ! was single hitting single + 2.0_pReal * (abs(rhoSgl(1:ns,2*c+3)) * abs(gdot(1:ns,2*c)) + abs(rhoSgl(1:ns,2*c+4)) * abs(gdot(1:ns,2*c-1))) & ! was single hitting immobile single or was immobile single hit by single @@ -2917,8 +2917,8 @@ if (lattice_structure(ph) == LATTICE_fcc_ID) & ! only fcc rhoDotAthermalAnnihilation(colinearSystem(s,instance),1:2) = - rhoDotAthermalAnnihilation(s,10) & * 0.25_pReal * sqrt(rhoForest(s)) * (dUpper(s,2) + dLower(s,2)) * edgeJogFactor(instance) - - + + !*** thermally activated annihilation of edge dipoles by climb rhoDotThermalAnnihilation = 0.0_pReal @@ -2942,7 +2942,7 @@ rhoDot = rhoDotFlux & + rhoDotMultiplication & + rhoDotSingle2DipoleGlide & + rhoDotAthermalAnnihilation & - + rhoDotThermalAnnihilation + + rhoDotThermalAnnihilation if (numerics_integrationMode == 1_pInt) then ! save rates for output if in central integration mode rhoDotFluxOutput(1:ns,1:8,1_pInt,ip,el) = rhoDotFlux(1:ns,1:8) @@ -2981,7 +2981,7 @@ endif if ( any(rhoSglOriginal(1:ns,1:4) + rhoDot(1:ns,1:4) * timestep < -aTolRho(instance)) & .or. any(rhoDipOriginal(1:ns,1:2) + rhoDot(1:ns,9:10) * timestep < -aTolRho(instance))) then #ifndef _OPENMP - if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt) then + if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt) then write(6,'(a,i5,a,i2)') '<< CONST >> evolution rate leads to negative density at el ',el,' ip ',ip write(6,'(a)') '<< CONST >> enforcing cutback !!!' endif @@ -2989,7 +2989,7 @@ if ( any(rhoSglOriginal(1:ns,1:4) + rhoDot(1:ns,1:4) * timestep < -aTolRho(in plasticState(p)%dotState = DAMASK_NaN return else - forall (s = 1:ns, t = 1_pInt:4_pInt) + forall (s = 1:ns, t = 1_pInt:4_pInt) plasticState(p)%dotState(iRhoU(s,t,instance),o) = rhoDot(s,t) plasticState(p)%dotState(iRhoB(s,t,instance),o) = rhoDot(s,t+4_pInt) endforall @@ -3037,10 +3037,10 @@ integer(pInt), intent(in) :: i, & e ! element index real(pReal), dimension(4,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & orientation ! crystal orientation in quaternions - + !* local variables integer(pInt) Nneighbors, & ! number of neighbors - n, & ! neighbor index + n, & ! neighbor index neighbor_e, & ! element index of my neighbor neighbor_i, & ! integration point index of my neighbor ph, & @@ -3054,15 +3054,15 @@ integer(pInt) Nneighbors, & real(pReal), dimension(4) :: absoluteMisorientation ! absolute misorientation (without symmetry) between me and my neighbor real(pReal), dimension(2,totalNslip(phase_plasticityInstance(material_phase(1,i,e))),& totalNslip(phase_plasticityInstance(material_phase(1,i,e))),& - FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,e))))) :: & + FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,e))))) :: & my_compatibility ! my_compatibility for current element and ip -real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(1,i,e)))) :: & +real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(1,i,e)))) :: & slipNormal, & slipDirection real(pReal) my_compatibilitySum, & thresholdValue, & nThresholdValues -logical, dimension(totalNslip(phase_plasticityInstance(material_phase(1,i,e)))) :: & +logical, dimension(totalNslip(phase_plasticityInstance(material_phase(1,i,e)))) :: & belowThreshold @@ -3087,24 +3087,24 @@ forall(s1 = 1_pInt:ns) & do n = 1_pInt,Nneighbors neighbor_e = mesh_ipNeighborhood(1,n,i,e) neighbor_i = mesh_ipNeighborhood(2,n,i,e) - - + + !* FREE SURFACE !* Set surface transmissivity to the value specified in the material.config - + if (neighbor_e <= 0_pInt .or. neighbor_i <= 0_pInt) then forall(s1 = 1_pInt:ns) & my_compatibility(1:2,s1,s1,n) = sqrt(surfaceTransmissivity(instance)) cycle endif - - + + !* PHASE BOUNDARY - !* If we encounter a different nonlocal "cpfem" phase at the neighbor, + !* If we encounter a different nonlocal "cpfem" phase at the neighbor, !* we consider this to be a real "physical" phase boundary, so completely incompatible. - !* If one of the two "CPFEM" phases has a local plasticity law, + !* If one of the two "CPFEM" phases has a local plasticity law, !* we do not consider this to be a phase boundary, so completely compatible. - + neighbor_phase = material_phase(1,neighbor_i,neighbor_e) if (neighbor_phase /= ph) then if (.not. phase_localPlasticity(neighbor_phase) .and. .not. phase_localPlasticity(ph)) then @@ -3114,7 +3114,7 @@ do n = 1_pInt,Nneighbors cycle endif - + !* GRAIN BOUNDARY ! !* fixed transmissivity for adjacent ips with different texture (only if explicitly given in material.config) @@ -3127,16 +3127,16 @@ do n = 1_pInt,Nneighbors endif cycle endif - + !* GRAIN BOUNDARY ? !* Compatibility defined by relative orientation of slip systems: !* The my_compatibility value is defined as the product of the slip normal projection and the slip direction projection. - !* Its sign is always positive for screws, for edges it has the same sign as the slip normal projection. - !* Since the sum for each slip system can easily exceed one (which would result in a transmissivity larger than one), + !* Its sign is always positive for screws, for edges it has the same sign as the slip normal projection. + !* Since the sum for each slip system can easily exceed one (which would result in a transmissivity larger than one), !* only values above or equal to a certain threshold value are considered. This threshold value is chosen, such that - !* the number of compatible slip systems is minimized with the sum of the original my_compatibility values exceeding one. - !* Finally the smallest my_compatibility value is decreased until the sum is exactly equal to one. + !* the number of compatible slip systems is minimized with the sum of the original my_compatibility values exceeding one. + !* Finally the smallest my_compatibility value is decreased until the sum is exactly equal to one. !* All values below the threshold are set to zero. else absoluteMisorientation = lattice_qDisorientation(orientation(1:4,1,i,e), & @@ -3148,7 +3148,7 @@ do n = 1_pInt,Nneighbors my_compatibility(2,s2,s1,n) = abs(math_mul3x3(slipNormal(1:3,s1), math_qRot(absoluteMisorientation, slipNormal(1:3,s2)))) & * abs(math_mul3x3(slipDirection(1:3,s1), math_qRot(absoluteMisorientation, slipDirection(1:3,s2)))) enddo - + my_compatibilitySum = 0.0_pReal belowThreshold = .true. do while (my_compatibilitySum < 1.0_pReal .and. any(belowThreshold(1:ns))) @@ -3284,7 +3284,7 @@ if (.not. phase_localPlasticity(ph)) then invFe = math_inv33(Fe(1:3,1:3,1_pInt,ip,el)) !* in case of periodic surfaces we have to find out how many periodic images in each direction we need - + do dir = 1_pInt,3_pInt maxCoord(dir) = maxval(mesh_node0(dir,:)) minCoord(dir) = minval(mesh_node0(dir,:)) @@ -3299,10 +3299,10 @@ if (.not. phase_localPlasticity(ph)) then endif enddo - - !* loop through all material points (also through their periodic images if present), + + !* loop through all material points (also through their periodic images if present), !* but only consider nonlocal neighbors within a certain cutoff radius R - + do neighbor_el = 1_pInt,mesh_NcpElems ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el))) neighbor_phase = material_phase(1_pInt,neighbor_ip,neighbor_el) @@ -3314,7 +3314,7 @@ if (.not. phase_localPlasticity(ph)) then neighbor_ns = totalNslip(neighbor_instance) neighbor_invFe = math_inv33(Fe(1:3,1:3,1,neighbor_ip,neighbor_el)) neighbor_ipVolumeSideLength = mesh_ipVolume(neighbor_ip,neighbor_el) ** (1.0_pReal/3.0_pReal) ! reference volume used here - + forall (s = 1_pInt:neighbor_ns, c = 1_pInt:2_pInt) neighbor_rhoExcess(c,1,s) = plasticState(np)%state(iRhoU(s,2*c-1,neighbor_instance),no) & ! positive mobiles - plasticState(np)%state(iRhoU(s,2*c,neighbor_instance),no) ! negative mobiles @@ -3326,43 +3326,43 @@ if (.not. phase_localPlasticity(ph)) then do deltaX = periodicImages(1,1),periodicImages(2,1) do deltaY = periodicImages(1,2),periodicImages(2,2) do deltaZ = periodicImages(1,3),periodicImages(2,3) - - + + !* regular case - + if (neighbor_el /= el .or. neighbor_ip /= ip & .or. deltaX /= 0_pInt .or. deltaY /= 0_pInt .or. deltaZ /= 0_pInt) then - + neighbor_coords = mesh_cellCenterCoordinates(neighbor_ip,neighbor_el) & + [real(deltaX,pReal), real(deltaY,pReal), real(deltaZ,pReal)] * meshSize connection = neighbor_coords - coords distance = sqrt(sum(connection * connection)) if (distance > cutoffRadius(instance)) cycle - + !* the segment length is the minimum of the third root of the control volume and the ip distance !* this ensures, that the central MP never sits on a neighbor dislocation segment - + connection_neighborLattice = math_mul33x3(neighbor_invFe, connection) segmentLength = min(neighbor_ipVolumeSideLength, distance) - + !* loop through all slip systems of the neighbor material point !* and add up the stress contributions from egde and screw excess on these slip systems (if significant) - + do s = 1_pInt,neighbor_ns if (all(abs(neighbor_rhoExcess(:,:,s)) < significantRho(instance))) cycle ! not significant - - + + !* map the connection vector from the lattice into the slip system frame - + connection_neighborSlip = math_mul33x3(lattice2slip(1:3,1:3,s,neighbor_instance), & connection_neighborLattice) - - + + !* edge contribution to stress sigma = 0.0_pReal - + x = connection_neighborSlip(1) y = connection_neighborSlip(2) z = connection_neighborSlip(3) @@ -3372,7 +3372,7 @@ if (.not. phase_localPlasticity(ph)) then do j = 1_pInt,2_pInt if (abs(neighbor_rhoExcess(1,j,s)) < significantRho(instance)) then - cycle + cycle elseif (j > 1_pInt) then x = connection_neighborSlip(1) & + sign(0.5_pReal * segmentLength, & @@ -3381,16 +3381,16 @@ if (.not. phase_localPlasticity(ph)) then xsquare = x * x endif - + flipSign = sign(1.0_pReal, -y) do side = 1_pInt,-1_pInt,-2_pInt lambda = real(side,pReal) * 0.5_pReal * segmentLength - y R = sqrt(xsquare + zsquare + lambda * lambda) Rsquare = R * R - Rcube = Rsquare * R + Rcube = Rsquare * R denominator = R * (R + flipSign * lambda) if (abs(denominator)<= tiny(0.0_pReal)) exit ipLoop - + sigma(1,1) = sigma(1,1) - real(side,pReal) & * flipSign * z / denominator & * (1.0_pReal + xsquare / Rsquare + xsquare / denominator) & @@ -3411,14 +3411,14 @@ if (.not. phase_localPlasticity(ph)) then sigma(2,3) = sigma(2,3) - real(side,pReal) & * (lattice_nu(ph) / R - zsquare / Rcube) * neighbor_rhoExcess(1,j,s) enddo - enddo - + enddo + !* screw contribution to stress - + x = connection_neighborSlip(1) ! have to restore this value, because position might have been adapted for edge deads before do j = 1_pInt,2_pInt if (abs(neighbor_rhoExcess(2,j,s)) < significantRho(instance)) then - cycle + cycle elseif (j > 1_pInt) then y = connection_neighborSlip(2) & + sign(0.5_pReal * segmentLength, & @@ -3432,10 +3432,10 @@ if (.not. phase_localPlasticity(ph)) then lambda = x + real(side,pReal) * 0.5_pReal * segmentLength R = sqrt(ysquare + zsquare + lambda * lambda) Rsquare = R * R - Rcube = Rsquare * R + Rcube = Rsquare * R denominator = R * (R + flipSign * lambda) if (abs(denominator)<= tiny(0.0_pReal)) exit ipLoop - + sigma(1,2) = sigma(1,2) - real(side,pReal) * flipSign * z & * (1.0_pReal - lattice_nu(ph)) / denominator & * neighbor_rhoExcess(2,j,s) @@ -3444,25 +3444,25 @@ if (.not. phase_localPlasticity(ph)) then * neighbor_rhoExcess(2,j,s) enddo enddo - + if (all(abs(sigma) < 1.0e-10_pReal)) cycle ! SIGMA IS NOT A REAL STRESS, THATS WHY WE NEED A REALLY SMALL VALUE HERE !* copy symmetric parts - + sigma(2,1) = sigma(1,2) sigma(3,1) = sigma(1,3) sigma(3,2) = sigma(2,3) - + !* scale stresses and map them into the neighbor material point's lattice configuration - + sigma = sigma * lattice_mu(neighbor_phase) * burgers(s,neighbor_instance) & / (4.0_pReal * pi * (1.0_pReal - lattice_nu(neighbor_phase))) & * mesh_ipVolume(neighbor_ip,neighbor_el) / segmentLength ! reference volume is used here (according to the segment length calculation) Tdislo_neighborLattice = Tdislo_neighborLattice & + math_mul33x33(math_transpose33(lattice2slip(1:3,1:3,s,neighbor_instance)), & math_mul33x33(sigma, lattice2slip(1:3,1:3,s,neighbor_instance))) - + enddo ! slip system loop @@ -3470,16 +3470,16 @@ if (.not. phase_localPlasticity(ph)) then !* only consider dead dislocations !* we assume that they all sit at a distance equal to half the third root of V !* in direction of the according slip direction - + else - + forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) & - rhoExcessDead(c,s) = plasticState(p)%state(iRhoB(s,2*c-1,instance),o) & ! positive deads (here we use symmetry: if this has negative sign it is - !treated as negative density at positive position instead of positive + rhoExcessDead(c,s) = plasticState(p)%state(iRhoB(s,2*c-1,instance),o) & ! positive deads (here we use symmetry: if this has negative sign it is + !treated as negative density at positive position instead of positive !density at negative position) - + plasticState(p)%state(iRhoB(s,2*c,instance),o) ! negative deads (here we use symmetry: if this has negative sign it is - !treated as positive density at positive position instead of negative + + plasticState(p)%state(iRhoB(s,2*c,instance),o) ! negative deads (here we use symmetry: if this has negative sign it is + !treated as positive density at positive position instead of negative !density at negative position) do s = 1_pInt,ns if (all(abs(rhoExcessDead(:,s)) < significantRho(instance))) cycle ! not significant @@ -3488,11 +3488,11 @@ if (.not. phase_localPlasticity(ph)) then * neighbor_ipVolumeSideLength * lattice_mu(ph) * burgers(s,instance) & / (sqrt(2.0_pReal) * pi * (1.0_pReal - lattice_nu(ph))) sigma(3,1) = sigma(1,3) - + Tdislo_neighborLattice = Tdislo_neighborLattice & + math_mul33x33(math_transpose33(lattice2slip(1:3,1:3,s,instance)), & math_mul33x33(sigma, lattice2slip(1:3,1:3,s,instance))) - + enddo ! slip system loop endif @@ -3502,7 +3502,7 @@ if (.not. phase_localPlasticity(ph)) then enddo ! deltaX loop - !* map the stress from the neighbor MP's lattice configuration into the deformed configuration + !* map the stress from the neighbor MP's lattice configuration into the deformed configuration !* and back into my lattice configuration neighborLattice2myLattice = math_mul33x33(invFe, Fe(1:3,1:3,1,neighbor_ip,neighbor_el)) @@ -3510,15 +3510,15 @@ if (.not. phase_localPlasticity(ph)) then + math_mul33x33(neighborLattice2myLattice, & math_mul33x33(Tdislo_neighborLattice, & math_transpose33(neighborLattice2myLattice))) - + enddo ipLoop enddo ! element loop - + endif end function plastic_nonlocal_dislocationstress - + !-------------------------------------------------------------------------------------------------- !> @brief return array of constitutive results !-------------------------------------------------------------------------------------------------- @@ -3553,11 +3553,11 @@ function plastic_nonlocal_postResults(Tstar_v,Fe,ip,el) integer(pInt), intent(in) :: & ip, & !< integration point el !< element - + real(pReal), dimension(plastic_nonlocal_sizePostResults(& phase_plasticityInstance(material_phase(1_pInt,ip,el)))) :: & plastic_nonlocal_postResults - + integer(pInt) :: & ph, & instance, & !< current instance of this plasticity @@ -3623,7 +3623,7 @@ tauBack = plasticState(ph)%State(iTauB(1:ns,instance),of) forall (t = 1_pInt:4_pInt) & gdot(1:ns,t) = rhoSgl(1:ns,t) * burgers(1:ns,instance) * v(1:ns,t) - + !* calculate limits for stable dipole height @@ -3641,7 +3641,7 @@ dUpper(1:ns,2) = lattice_mu(ph) * burgers(1:ns,instance) & forall (c = 1_pInt:2_pInt) where(sqrt(rhoSgl(1:ns,2*c-1)+rhoSgl(1:ns,2*c)+& abs(rhoSgl(1:ns,2*c+3))+abs(rhoSgl(1:ns,2*c+4))+rhoDip(1:ns,c)) >= tiny(0.0_pReal)) & - dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) & + dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) & + abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), & dUpper(1:ns,c)) end forall @@ -3664,95 +3664,95 @@ outputsLoop: do o = 1_pInt,plastic_nonlocal_Noutput(instance) case (rho_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(abs(rhoSgl),2) + sum(rhoDip,2) cs = cs + ns - + case (rho_sgl_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(abs(rhoSgl),2) cs = cs + ns - + case (rho_sgl_mobile_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(abs(rhoSgl(1:ns,1:4)),2) cs = cs + ns - + case (rho_sgl_immobile_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoSgl(1:ns,5:8),2) cs = cs + ns - + case (rho_dip_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDip,2) cs = cs + ns - + case (rho_edge_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(abs(rhoSgl(1:ns,[1,2,5,6])),2) + rhoDip(1:ns,1) cs = cs + ns - + case (rho_sgl_edge_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(abs(rhoSgl(1:ns,[1,2,5,6])),2) cs = cs + ns - + case (rho_sgl_edge_mobile_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoSgl(1:ns,1:2),2) cs = cs + ns - + case (rho_sgl_edge_immobile_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoSgl(1:ns,5:6),2) cs = cs + ns - + case (rho_sgl_edge_pos_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,1) + abs(rhoSgl(1:ns,5)) cs = cs + ns - + case (rho_sgl_edge_pos_mobile_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,1) cs = cs + ns - + case (rho_sgl_edge_pos_immobile_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,5) cs = cs + ns - + case (rho_sgl_edge_neg_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,2) + abs(rhoSgl(1:ns,6)) cs = cs + ns - + case (rho_sgl_edge_neg_mobile_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,2) cs = cs + ns - + case (rho_sgl_edge_neg_immobile_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,6) cs = cs + ns - + case (rho_dip_edge_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDip(1:ns,1) cs = cs + ns - + case (rho_screw_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(abs(rhoSgl(1:ns,[3,4,7,8])),2) + rhoDip(1:ns,2) cs = cs + ns - + case (rho_sgl_screw_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(abs(rhoSgl(1:ns,[3,4,7,8])),2) cs = cs + ns - + case (rho_sgl_screw_mobile_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoSgl(1:ns,3:4),2) cs = cs + ns - + case (rho_sgl_screw_immobile_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoSgl(1:ns,7:8),2) cs = cs + ns - + case (rho_sgl_screw_pos_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,3) + abs(rhoSgl(1:ns,7)) cs = cs + ns - + case (rho_sgl_screw_pos_mobile_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,3) cs = cs + ns - + case (rho_sgl_screw_pos_immobile_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,7) cs = cs + ns - + case (rho_sgl_screw_neg_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,4) + abs(rhoSgl(1:ns,8)) cs = cs + ns @@ -3768,82 +3768,82 @@ outputsLoop: do o = 1_pInt,plastic_nonlocal_Noutput(instance) case (rho_dip_screw_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDip(1:ns,2) cs = cs + ns - + case (excess_rho_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = (rhoSgl(1:ns,1) + abs(rhoSgl(1:ns,5))) & - (rhoSgl(1:ns,2) + abs(rhoSgl(1:ns,6))) & + (rhoSgl(1:ns,3) + abs(rhoSgl(1:ns,7))) & - (rhoSgl(1:ns,4) + abs(rhoSgl(1:ns,8))) cs = cs + ns - + case (excess_rho_edge_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = (rhoSgl(1:ns,1) + abs(rhoSgl(1:ns,5))) & - (rhoSgl(1:ns,2) + abs(rhoSgl(1:ns,6))) cs = cs + ns - + case (excess_rho_screw_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = (rhoSgl(1:ns,3) + abs(rhoSgl(1:ns,7))) & - (rhoSgl(1:ns,4) + abs(rhoSgl(1:ns,8))) cs = cs + ns - + case (rho_forest_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoForest cs = cs + ns - + case (delta_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = 1.0_pReal / sqrt(sum(abs(rhoSgl),2) + sum(rhoDip,2)) cs = cs + ns - + case (delta_sgl_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = 1.0_pReal / sqrt(sum(abs(rhoSgl),2)) cs = cs + ns - + case (delta_dip_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = 1.0_pReal / sqrt(sum(rhoDip,2)) cs = cs + ns - + case (shearrate_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(gdot,2) cs = cs + ns - + case (resolvedstress_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = tau cs = cs + ns - + case (resolvedstress_back_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = tauBack cs = cs + ns - + case (resolvedstress_external_ID) - do s = 1_pInt,ns + do s = 1_pInt,ns sLattice = slipSystemLattice(s,instance) plastic_nonlocal_postResults(cs+s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,ph)) enddo cs = cs + ns - + case (resistance_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = tauThreshold cs = cs + ns - + case (rho_dot_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotSgl(1:ns,1:4),2) & + sum(rhoDotSgl(1:ns,5:8)*sign(1.0_pReal,rhoSgl(1:ns,5:8)),2) & + sum(rhoDotDip,2) cs = cs + ns - + case (rho_dot_sgl_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotSgl(1:ns,1:4),2) & + sum(rhoDotSgl(1:ns,5:8)*sign(1.0_pReal,rhoSgl(1:ns,5:8)),2) cs = cs + ns - + case (rho_dot_sgl_mobile_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotSgl(1:ns,1:4),2) cs = cs + ns - + case (rho_dot_dip_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotDip,2) cs = cs + ns - + case (rho_dot_gen_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotMultiplicationOutput(1:ns,1,1_pInt,ip,el) & + rhoDotMultiplicationOutput(1:ns,2,1_pInt,ip,el) @@ -3856,157 +3856,157 @@ outputsLoop: do o = 1_pInt,plastic_nonlocal_Noutput(instance) case (rho_dot_gen_screw_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotMultiplicationOutput(1:ns,2,1_pInt,ip,el) cs = cs + ns - + case (rho_dot_sgl2dip_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotSingle2DipoleGlideOutput(1:ns,1,1_pInt,ip,el) & + rhoDotSingle2DipoleGlideOutput(1:ns,2,1_pInt,ip,el) cs = cs + ns - + case (rho_dot_sgl2dip_edge_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotSingle2DipoleGlideOutput(1:ns,1,1_pInt,ip,el) cs = cs + ns - + case (rho_dot_sgl2dip_screw_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotSingle2DipoleGlideOutput(1:ns,2,1_pInt,ip,el) cs = cs + ns - + case (rho_dot_ann_ath_ID) - plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotAthermalAnnihilationOutput(1:ns,1,1_pInt,ip,el) & + plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotAthermalAnnihilationOutput(1:ns,1,1_pInt,ip,el) & + rhoDotAthermalAnnihilationOutput(1:ns,2,1_pInt,ip,el) cs = cs + ns - - case (rho_dot_ann_the_ID) - plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotThermalAnnihilationOutput(1:ns,1,1_pInt,ip,el) & + + case (rho_dot_ann_the_ID) + plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotThermalAnnihilationOutput(1:ns,1,1_pInt,ip,el) & + rhoDotThermalAnnihilationOutput(1:ns,2,1_pInt,ip,el) cs = cs + ns - case (rho_dot_ann_the_edge_ID) - plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotThermalAnnihilationOutput(1:ns,1,1_pInt,ip,el) + case (rho_dot_ann_the_edge_ID) + plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotThermalAnnihilationOutput(1:ns,1,1_pInt,ip,el) cs = cs + ns - case (rho_dot_ann_the_screw_ID) + case (rho_dot_ann_the_screw_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotThermalAnnihilationOutput(1:ns,2,1_pInt,ip,el) cs = cs + ns - case (rho_dot_edgejogs_ID) + case (rho_dot_edgejogs_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotEdgeJogsOutput(1:ns,1_pInt,ip,el) cs = cs + ns case (rho_dot_flux_mobile_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotFluxOutput(1:ns,1:4,1_pInt,ip,el),2) cs = cs + ns - + case (rho_dot_flux_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotFluxOutput(1:ns,1:4,1_pInt,ip,el),2) & + sum(rhoDotFluxOutput(1:ns,5:8,1_pInt,ip,el)*sign(1.0_pReal,rhoSgl(1:ns,5:8)),2) cs = cs + ns - + case (rho_dot_flux_edge_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotFluxOutput(1:ns,1:2,1_pInt,ip,el),2) & + sum(rhoDotFluxOutput(1:ns,5:6,1_pInt,ip,el)*sign(1.0_pReal,rhoSgl(1:ns,5:6)),2) cs = cs + ns - + case (rho_dot_flux_screw_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotFluxOutput(1:ns,3:4,1_pInt,ip,el),2) & + sum(rhoDotFluxOutput(1:ns,7:8,1_pInt,ip,el)*sign(1.0_pReal,rhoSgl(1:ns,7:8)),2) cs = cs + ns - + case (velocity_edge_pos_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = v(1:ns,1) cs = cs + ns - + case (velocity_edge_neg_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = v(1:ns,2) cs = cs + ns - + case (velocity_screw_pos_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = v(1:ns,3) cs = cs + ns - + case (velocity_screw_neg_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = v(1:ns,4) cs = cs + ns - + case (slipdirectionx_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = m_currentconf(1,1:ns,1) cs = cs + ns - + case (slipdirectiony_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = m_currentconf(2,1:ns,1) cs = cs + ns - + case (slipdirectionz_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = m_currentconf(3,1:ns,1) cs = cs + ns - + case (slipnormalx_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = n_currentconf(1,1:ns) cs = cs + ns - + case (slipnormaly_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = n_currentconf(2,1:ns) cs = cs + ns - + case (slipnormalz_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = n_currentconf(3,1:ns) cs = cs + ns - + case (fluxdensity_edge_posx_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,1) * v(1:ns,1) * m_currentconf(1,1:ns,1) cs = cs + ns - + case (fluxdensity_edge_posy_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,1) * v(1:ns,1) * m_currentconf(2,1:ns,1) cs = cs + ns - + case (fluxdensity_edge_posz_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,1) * v(1:ns,1) * m_currentconf(3,1:ns,1) cs = cs + ns - + case (fluxdensity_edge_negx_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = - rhoSgl(1:ns,2) * v(1:ns,2) * m_currentconf(1,1:ns,1) cs = cs + ns - + case (fluxdensity_edge_negy_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = - rhoSgl(1:ns,2) * v(1:ns,2) * m_currentconf(2,1:ns,1) cs = cs + ns - + case (fluxdensity_edge_negz_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = - rhoSgl(1:ns,2) * v(1:ns,2) * m_currentconf(3,1:ns,1) cs = cs + ns - + case (fluxdensity_screw_posx_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,3) * v(1:ns,3) * m_currentconf(1,1:ns,2) cs = cs + ns - + case (fluxdensity_screw_posy_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,3) * v(1:ns,3) * m_currentconf(2,1:ns,2) cs = cs + ns - + case (fluxdensity_screw_posz_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,3) * v(1:ns,3) * m_currentconf(3,1:ns,2) cs = cs + ns - + case (fluxdensity_screw_negx_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = - rhoSgl(1:ns,4) * v(1:ns,4) * m_currentconf(1,1:ns,2) cs = cs + ns - + case (fluxdensity_screw_negy_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = - rhoSgl(1:ns,4) * v(1:ns,4) * m_currentconf(2,1:ns,2) cs = cs + ns - + case (fluxdensity_screw_negz_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = - rhoSgl(1:ns,4) * v(1:ns,4) * m_currentconf(3,1:ns,2) cs = cs + ns - + case (maximumdipoleheight_edge_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = dUpper(1:ns,1) cs = cs + ns - + case (maximumdipoleheight_screw_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = dUpper(1:ns,2) cs = cs + ns - + case(dislocationstress_ID) sigma = plastic_nonlocal_dislocationstress(Fe, ip, el) plastic_nonlocal_postResults(cs+1_pInt) = sigma(1,1) @@ -4016,11 +4016,11 @@ outputsLoop: do o = 1_pInt,plastic_nonlocal_Noutput(instance) plastic_nonlocal_postResults(cs+5_pInt) = sigma(2,3) plastic_nonlocal_postResults(cs+6_pInt) = sigma(3,1) cs = cs + 6_pInt - + case(accumulatedshear_ID) plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = plasticState(ph)%state(iGamma(1:ns,instance),of) cs = cs + ns - + end select enddo outputsLoop diff --git a/src/plastic_phenoplus.f90 b/src/plastic_phenoplus.f90 index 0887da239..38561af60 100644 --- a/src/plastic_phenoplus.f90 +++ b/src/plastic_phenoplus.f90 @@ -557,7 +557,7 @@ subroutine plastic_phenoplus_init(fileUnit) ! allocate state arrays sizeState = plastic_phenoplus_totalNslip(instance) & ! s_slip + plastic_phenoplus_totalNtwin(instance) & ! s_twin - + 2_pInt & ! sum(gamma) + sum(f) + + 2_pInt & ! sum(gamma) + sum(twinVolFrac) + plastic_phenoplus_totalNslip(instance) & ! accshear_slip + plastic_phenoplus_totalNtwin(instance) & ! accshear_twin + plastic_phenoplus_totalNslip(instance) ! kappa @@ -568,7 +568,7 @@ subroutine plastic_phenoplus_init(fileUnit) ! memory leak issue. sizeDotState = plastic_phenoplus_totalNslip(instance) & ! s_slip + plastic_phenoplus_totalNtwin(instance) & ! s_twin - + 2_pInt & ! sum(gamma) + sum(f) + + 2_pInt & ! sum(gamma) + sum(twinVolFrac) + plastic_phenoplus_totalNslip(instance) & ! accshear_slip + plastic_phenoplus_totalNtwin(instance) ! accshear_twin @@ -739,183 +739,256 @@ end subroutine plastic_phenoplus_aTolState !-------------------------------------------------------------------------------------------------- !> @brief calculate push-up factors (kappa) for each voxel based on its neighbors !-------------------------------------------------------------------------------------------------- -subroutine plastic_phenoplus_microstructure(orientation,ipc,ip,el) - use math, only: pi, & - math_mul33x33, & - math_mul3x3, & - math_transpose33, & - math_qDot, & - math_qRot, & - indeg +subroutine plastic_phenoplus_microstructure(orientation,ipc,ip,el,F0,Fe,Fp,Tstar_v) + use math, only: pi, & + math_identity2nd, & + math_mul33x33, & + math_mul33xx33, & + math_mul3x3, & + math_transpose33, & + math_qDot, & + math_qRot, & + indeg - use mesh, only: mesh_element, & - FE_NipNeighbors, & - FE_geomtype, & - FE_celltype, & - mesh_maxNips, & - mesh_NcpElems, & - mesh_ipNeighborhood + use mesh, only: mesh_element, & + FE_NipNeighbors, & + FE_geomtype, & + FE_celltype, & + mesh_maxNips, & + mesh_NcpElems, & + mesh_ipNeighborhood - use material, only: material_phase, & - material_texture, & - phase_plasticityInstance, & - phaseAt, phasememberAt, & - homogenization_maxNgrains, & - plasticState + use material, only: material_phase, & + material_texture, & + phase_plasticityInstance, & + phaseAt, phasememberAt, & + homogenization_maxNgrains, & + plasticState - use lattice, only: lattice_sn, & - lattice_sd, & - lattice_qDisorientation + use lattice, only: lattice_Sslip_v, & + lattice_maxNslipFamily, & + lattice_NslipSystem, & + lattice_NslipSystem, & + lattice_sn, & + lattice_sd, & + lattice_qDisorientation - !***input variables - implicit none - integer(pInt), intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - real(pReal), dimension(4,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & - orientation ! crystal orientation in quaternions + !***input variables + implicit none + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el + real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + F0, & !< deformation gradient from last increment + Fe, & !< elastic deformation gradient + Fp !< elastic deformation gradient !< element + real(pReal), dimension(4,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + orientation !< crystal orientation in quaternions + real(pReal), dimension(6,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + Tstar_v !< for calculation of gdot - !***local variables - integer(pInt) instance, & !my instance of this plasticity - ph, & !my phase - of, & !my spatial position in memory (offset) - textureID, & !my texture - Nneighbors, & !number of neighbors (<= 6) - vld_Nneighbors, & !number of my valid neighbors - n, & !neighbor index (for iterating through all neighbors) - ns, & !number of slip system - nt, & !number of twin system - me_slip, & !my slip system index - neighbor_el, & !element number of neighboring material point - neighbor_ip, & !integration point of neighboring material point - neighbor_n, & !I have no idea what is this - neighbor_of, & !spatial position in memory for this neighbor (offset) - neighbor_ph, & !neighbor's phase - neighbor_tex, & !neighbor's texture ID - ne_slip_ac, & !loop to find neighbor shear - ne_slip, & !slip system index for neighbor - index_kappa, & !index of pushup factors in plasticState - offset_acshear_slip, & !offset in PlasticState for the accumulative shear - j !quickly loop through slip families + !***local variables + integer(pInt) instance, & !my instance of this plasticity + ph, & !my phase + of, & !my spatial position in memory (offset) + textureID, & !my texture + index_myFamily, & + Nneighbors, & !number of neighbors (<= 6) + vld_Nneighbors, & !number of my valid neighbors + n, & !neighbor index (for iterating through all neighbors) + n_calcTaylor, & ! + n_phasecheck, & ! + ns, & !number of slip system + nt, & !number of twin system + me_slip, & !my slip system index + neighbor_el, & !element number of neighboring material point + neighbor_ip, & !integration point of neighboring material point + neighbor_ipc, & !I have no idea what is this + neighbor_of, & !spatial position in memory for this neighbor (offset) + neighbor_ph, & !neighbor's phase + neighbor_instance, & !neighbor's instance of this plasticity + neighbor_tex, & !neighbor's texture ID + ne_slip, & !slip system index for neighbor + index_kappa, & !index of pushup factors in plasticState + j, & !quickly loop through slip families + f,i,& !loop counter for me + f_ne, i_ne !loop counter for neighbor - real(pReal) kappa_max, & ! - tmp_myshear_slip, & !temp storage for accumulative shear for me - mprime_cut, & !m' cutoff to consider neighboring effect - avg_acshear_ne, & !the average accumulative shear from my neighbor - tmp_mprime, & !temp holder for m' value - tmp_acshear !temp holder for accumulative shear for m' + real(pReal) mprime_cut, & !m' cutoff to consider neighboring effect + dtaylor_cut, & !threshold for determine high contrast interface using Taylor factor + tau_slip, & !the average accumulative shear from my neighbor + taylor_me, & !Taylor factor for me + taylor_ne, & !Taylor factor for my current neighbor + d_vonstrain, & !von Mises delta strain (temp container) + sum_gdot !total shear rate for given material point + real(pReal), dimension(3,3) :: & + F0_me, & !my deformation gradient from last converged increment + Fe_me, & !my elastic deformation gradient + Fp_me, & !my plastic deformation gradient + dF_me, & !my deformation gradient change (delta) + dE_me, & !my Green Lagrangian strain tensor (delta) + F0_ne, & ! + Fe_ne, & !elastic deformation gradient of my current neighbor + Fp_ne, & !plastic deformation gradient of my current neighbor + dF_ne, & !deformation gradient of my current neighbor + dE_ne !delta Green Lagrangian strain tensor - real(pReal), dimension(plastic_phenoplus_totalNslip(phase_plasticityInstance(material_phase(1,ip,el)))) :: & - m_primes, & !m' between me_alpha(one) and neighbor beta(all) - me_acshear, & !temp storage for ac_shear of one particular system for me - ne_acshear !temp storage for ac_shear of one particular system for one of my neighbor + real(pReal), dimension(plastic_phenoplus_totalNslip(phase_plasticityInstance(material_phase(1,ip,el)))) :: & + m_primes !m' between me_alpha(one) and neighbor beta(all) - real(pReal), dimension(3,plastic_phenoplus_totalNslip(phase_plasticityInstance(material_phase(1,ip,el)))) :: & - slipNormal, & - slipDirect + real(pReal), dimension(3,plastic_phenoplus_totalNslip(phase_plasticityInstance(material_phase(1,ip,el)))) :: & + slipNormal, & + slipDirect - real(pReal), dimension(4) :: my_orientation, & !store my orientation - neighbor_orientation, & !store my neighbor orientation - absMisorientation + real(pReal), dimension(4) :: & + my_orientation, & !store my orientation + neighbor_orientation, & !store my neighbor orientation + absMisorientation - real(pReal), dimension(FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el))))) :: & - ne_mprimes !m' between each neighbor + real(pReal), dimension(FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el))))) :: & + ne_mprimes, & !m' between each neighbor + d_taylors !store (taylor_ne-taylor_me) for each neighbor - !***Get my properties - Nneighbors = FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el)))) - ph = phaseAt(ipc,ip,el) !get my phase - of = phasememberAt(ipc,ip,el) !get my spatial location offset in memory - textureID = material_texture(1,ip,el) !get my texture ID - instance = phase_plasticityInstance(ph) !get my instance based on phase ID - ns = plastic_phenoplus_totalNslip(instance) - nt = plastic_phenoplus_totalNtwin(instance) - offset_acshear_slip = ns + nt + 2_pInt - index_kappa = ns + nt + 2_pInt + ns + nt !location of kappa in plasticState - mprime_cut = 0.7_pReal !set by Dr.Bieler + !***Get my properties + Nneighbors = FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el)))) + ph = phaseAt(ipc,ip,el) !get my phase + of = phasememberAt(ipc,ip,el) !get my spatial location offset in memory + textureID = material_texture(1,ip,el) !get my texture ID + instance = phase_plasticityInstance(ph) !get my instance based on phase ID + ns = plastic_phenoplus_totalNslip(instance) + nt = plastic_phenoplus_totalNtwin(instance) + index_kappa = ns + nt + 2_pInt + ns + nt !location of kappa in plasticState - !***gather my accumulative shear from palsticState - FINDMYSHEAR: do j = 1_pInt,ns - me_acshear(j) = plasticState(ph)%state(offset_acshear_slip+j, of) - enddo FINDMYSHEAR + !***init calculation for given voxel + mprime_cut = 0.7_pReal !set by Dr.Bieler + dtaylor_cut = 1.0_pReal !set by Chen, quick test only - !***gather my orientation and slip systems - my_orientation = orientation(1:4, ipc, ip, el) - slipNormal(1:3, 1:ns) = lattice_sn(1:3, 1:ns, ph) - slipDirect(1:3, 1:ns) = lattice_sd(1:3, 1:ns, ph) - kappa_max = plastic_phenoplus_kappa_max(instance) !maximum pushups allowed (READIN) + !***gather my orientation, F and slip systems + my_orientation = orientation(1:4, ipc, ip, el) + F0_me = F0(1:3, 1:3, ipc, ip, el) + Fe_me = Fe(1:3, 1:3, ipc, ip, el) + Fp_me = Fp(1:3, 1:3, ipc, ip, el) + slipNormal(1:3, 1:ns) = lattice_sn(1:3, 1:ns, ph) + slipDirect(1:3, 1:ns) = lattice_sd(1:3, 1:ns, ph) - !***calculate kappa between me and all my neighbors - LOOPMYSLIP: DO me_slip=1_pInt,ns - vld_Nneighbors = Nneighbors - tmp_myshear_slip = me_acshear(me_slip) - tmp_mprime = 0.0_pReal !highest m' from all neighbors - tmp_acshear = 0.0_pReal !accumulative shear from highest m' + !***check if all my neighbors have the same phase as me + vld_Nneighbors = 0 + PHASECHECK: DO n_phasecheck = 1_pInt, Nneighbors + !******for each of my neighbor + neighbor_el = mesh_ipNeighborhood( 1, n_phasecheck, ip, el ) + neighbor_ip = mesh_ipNeighborhood( 2, n_phasecheck, ip, el ) + neighbor_ipc = 1 + neighbor_of = phasememberAt( neighbor_ipc, neighbor_ip, neighbor_el ) + neighbor_ph = phaseAt( neighbor_ipc, neighbor_ip, neighbor_el ) + IF (neighbor_ph == ph) THEN + vld_Nneighbors = vld_Nneighbors + 1_pInt + ENDIF + ENDDO PHASECHECK - !***go through my neighbors to find highest m' - LOOPNEIGHBORS: DO n=1_pInt,Nneighbors - neighbor_el = mesh_ipNeighborhood(1,n,ip,el) - neighbor_ip = mesh_ipNeighborhood(2,n,ip,el) - neighbor_n = 1 !It is ipc - neighbor_of = phasememberAt( neighbor_n, neighbor_ip, neighbor_el) - neighbor_ph = phaseAt( neighbor_n, neighbor_ip, neighbor_el) - neighbor_tex = material_texture(1,neighbor_ip,neighbor_el) - neighbor_orientation = orientation(1:4, neighbor_n, neighbor_ip, neighbor_el) !ipc is always 1. - absMisorientation = lattice_qDisorientation(my_orientation, & - neighbor_orientation, & - 0_pInt) !no need for explicit calculation of symmetry + !***initialize kappa with 1.0 (assume no push-up) + plasticState(ph)%state(index_kappa+1_pInt:index_kappa+ns, of) = 1.0_pReal - !***find the accumulative shear for this neighbor - LOOPFINDNEISHEAR: DO ne_slip_ac=1_pInt, ns - ne_acshear(ne_slip_ac) = plasticState(ph)%state(offset_acshear_slip+ne_slip_ac, & - neighbor_of) - ENDDO LOOPFINDNEISHEAR + !***only calculate kappa for those inside the main phase + IF (vld_Nneighbors == Nneighbors) THEN + !******calculate Taylor factor for me + dF_me = math_mul33x33(Fe_me,Fp_me) - F0_me + dE_me = 0.5*(math_mul33x33(math_transpose33(dF_me), dF_me) - math_identity2nd(3)) !dE = 0.5(dF^tdF-I) + d_vonstrain = SQRT(2.0_pReal/3.0_pReal * math_mul33xx33(dE_me, dE_me)) + sum_gdot = 0.0_pReal + !go through my slip system to find the sum of gamma_dot + j = 0_pInt + slipFamilies: DO f = 1_pInt,lattice_maxNslipFamily + index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) !at which index starts my family + slipSystems: DO i = 1_pInt,plastic_phenoplus_Nslip(f,instance) + j = j+1_pInt + tau_slip = dot_product(Tstar_v(1:6, ipc, ip, el),lattice_Sslip_v(1:6,1,index_myFamily+i,ph)) + sum_gdot = sum_gdot + & + plastic_phenoplus_gdot0_slip(instance)* & + ((abs(tau_slip)/(plasticState(ph)%state(j,of))) & + **plastic_phenoplus_n_slip(instance))*sign(1.0_pReal,tau_slip) + ENDDO slipSystems + ENDDO slipFamilies + taylor_me = d_vonstrain/sum_gdot - !***calculate the average accumulative shear and use it as cutoff - avg_acshear_ne = SUM(ne_acshear)/ns + !***calculate delta_M (Taylor factor) between each neighbor and me + LOOPCALCTAYLOR: DO n_calcTaylor=1_pInt, Nneighbors + !******for each of my neighbor + neighbor_el = mesh_ipNeighborhood( 1, n_calcTaylor, ip, el ) + neighbor_ip = mesh_ipNeighborhood( 2, n_calcTaylor, ip, el ) + neighbor_ipc = 1 !It is ipc + neighbor_of = phasememberAt( neighbor_ipc, neighbor_ip, neighbor_el ) + neighbor_ph = phaseAt( neighbor_ipc, neighbor_ip, neighbor_el ) + neighbor_instance = phase_plasticityInstance( neighbor_ph ) + neighbor_tex = material_texture( 1,neighbor_ip, neighbor_el ) + neighbor_orientation = orientation( 1:4, neighbor_ipc, neighbor_ip, neighbor_el ) !ipc is always 1. + Fe_ne = Fe( 1:3, 1:3, neighbor_ipc, neighbor_ip, neighbor_el ) + Fp_ne = Fp( 1:3, 1:3, neighbor_ipc, neighbor_ip, neighbor_el ) + F0_ne = F0( 1:3, 1:3, neighbor_ipc, neighbor_ip, neighbor_el ) + !******calculate the Taylor factor + dF_ne = math_mul33x33(Fe_ne, Fp_ne) - F0_ne + dE_ne = 0.5*(math_mul33x33(math_transpose33(dF_ne), dF_ne) - math_identity2nd(3)) !dE = 0.5(dF^tdF-I) + d_vonstrain = SQRT(2.0_pReal/3.0_pReal * math_mul33xx33(dE_ne, dE_ne)) + sum_gdot = 0.0_pReal + !go through my neighbor slip system to calculate sum_gdot + j = 0_pInt + slipFamiliesNeighbor: DO f_ne = 1_pInt,lattice_maxNslipFamily + index_myFamily = sum(lattice_NslipSystem(1:f_ne-1_pInt,neighbor_ph)) ! at which index starts my family + slipSystemsNeighbor: DO i_ne = 1_pInt,plastic_phenoplus_Nslip(f_ne,neighbor_instance) + j = j+1_pInt + tau_slip = dot_product(Tstar_v(1:6, neighbor_ipc, neighbor_ip, neighbor_el), & + lattice_Sslip_v(1:6,1,index_myFamily+i_ne,neighbor_ph)) + sum_gdot = sum_gdot & + +plastic_phenoplus_gdot0_slip(neighbor_instance) & + *((abs(tau_slip)/(plasticState(neighbor_ph)%state(j,neighbor_of))) & + **plastic_phenoplus_n_slip(neighbor_instance))*sign(1.0_pReal,tau_slip) + ENDDO slipSystemsNeighbor + ENDDO slipFamiliesNeighbor + taylor_ne = d_vonstrain / sum_gdot + !******calculate Taylor difference + d_taylors(n_calcTaylor) = taylor_ne - taylor_me + ENDDO LOOPCALCTAYLOR - !*** - IF (ph==neighbor_ph) THEN - !***walk through all the - LOOPNEIGHBORSLIP: DO ne_slip=1_pInt,ns - !***only consider slip system that is active (above average accumulative shear) - IF (ne_acshear(ne_slip) > avg_acshear_ne) THEN - m_primes(ne_slip) = abs(math_mul3x3(slipNormal(1:3,me_slip), & - math_qRot(absMisorientation, slipNormal(1:3,ne_slip)))) & - *abs(math_mul3x3(slipDirect(1:3,me_slip), & - math_qRot(absMisorientation, slipDirect(1:3,ne_slip)))) - !***find the highest m' and corresponding accumulative shear - IF (m_primes(ne_slip) > tmp_mprime) THEN - tmp_mprime = m_primes(ne_slip) - tmp_acshear = ne_acshear(ne_slip) - ENDIF - ENDIF - ENDDO LOOPNEIGHBORSLIP + !***Only perform necessary calculation if high contrast interface is detected + IF (maxval(d_taylors) > dtaylor_cut) THEN + !*****calculate kappa per slip system base + LOOPMYSLIP: DO me_slip = 1_pInt, ns + ne_mprimes = 0.0_pReal !initialize max m' to 0 for all neighbors + LOOPMYNEIGHBORS: DO n=1_pInt, Nneighbors + !*******only consider neighbor at the high contrast interface + IF (d_taylors(n) > dtaylor_cut) THEN + neighbor_el = mesh_ipNeighborhood( 1, n_calcTaylor, ip, el ) + neighbor_ip = mesh_ipNeighborhood( 2, n_calcTaylor, ip, el ) + neighbor_ipc = 1 !It is ipc + neighbor_of = phasememberAt( neighbor_ipc, neighbor_ip, neighbor_el ) + neighbor_ph = phaseAt( neighbor_ipc, neighbor_ip, neighbor_el ) + neighbor_instance = phase_plasticityInstance( neighbor_ph ) + neighbor_tex = material_texture( 1,neighbor_ip, neighbor_el ) + neighbor_orientation = orientation( 1:4, neighbor_ipc, neighbor_ip, neighbor_el ) !ipc is always 1. + absMisorientation = lattice_qDisorientation( my_orientation, & + neighbor_orientation, & + 0_pInt ) !no need for explicit calculation of symmetry + !*********go through neighbor slip system to calculate m' + LOOPNEIGHBORSLIP: DO ne_slip=1_pInt,ns + m_primes(ne_slip) = abs(math_mul3x3(slipNormal(1:3,me_slip), & + math_qRot(absMisorientation, slipNormal(1:3,ne_slip)))) & + *abs(math_mul3x3(slipDirect(1:3,me_slip), & + math_qRot(absMisorientation, slipDirect(1:3,ne_slip)))) + ENDDO LOOPNEIGHBORSLIP + ne_mprimes(n) = maxval(m_primes) + ENDIF + ENDDO LOOPMYNEIGHBORS - ELSE - ne_mprimes(n) = 0.0_pReal - vld_Nneighbors = vld_Nneighbors - 1_pInt - ENDIF + plasticState(ph)%state(index_kappa+me_slip, of) = & + 1.03_pReal + 0.03_pReal*ERF(4.0_pReal * maxval(ne_mprimes) - 4.0_pReal) - ENDDO LOOPNEIGHBORS + ENDDO LOOPMYSLIP + + ENDIF - !***check if this element close to rim - IF (vld_Nneighbors < Nneighbors) THEN - !***rim voxel, no modification allowed - plasticState(ph)%state(index_kappa+me_slip, of) = 1.0_pReal - ELSE - !***patch voxel, started to calculate push up factor for gamma_dot - IF ((tmp_mprime > mprime_cut) .AND. (tmp_acshear > tmp_myshear_slip)) THEN - plasticState(ph)%state(index_kappa+me_slip, of) = 1.0_pReal / tmp_mprime - ELSE - !***minimum damping factor is 0.5 - plasticState(ph)%state(index_kappa+me_slip, of) = 0.5_pReal + tmp_mprime * 0.5_pReal - ENDIF ENDIF - ENDDO LOOPMYSLIP - end subroutine plastic_phenoplus_microstructure