diff --git a/PRIVATE b/PRIVATE index 8fec909d1..174ecac2d 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 8fec909d1931b092b223b0560dd30c3339c6e5a7 +Subproject commit 174ecac2d3ab7596bdb60184d6bb9e1a52cb7378 diff --git a/examples/config/numerics.yaml b/examples/config/numerics.yaml index 72fbe82a2..4363b97cc 100644 --- a/examples/config/numerics.yaml +++ b/examples/config/numerics.yaml @@ -79,6 +79,5 @@ commercialFEM: unitlength: 1 # physical length of one computational length unit generic: - charLength: 1.0 # characteristic length scale for gradient problems. random_seed: 0 # fixed seeding for pseudo-random number generator, Default 0: use random seed. residualStiffness: 1.0e-6 # non-zero residual damage. diff --git a/examples/config/phase/damage/anisobrittle_cubic.yaml b/examples/config/phase/damage/anisobrittle_cubic.yaml index 5320a4aa9..41f64efcf 100644 --- a/examples/config/phase/damage/anisobrittle_cubic.yaml +++ b/examples/config/phase/damage/anisobrittle_cubic.yaml @@ -7,5 +7,5 @@ q: 20 output: [f_phi] -K_11: 1.0 +D_11: 1.0 mu: 0.001 diff --git a/examples/config/phase/damage/isobrittle_generic.yaml b/examples/config/phase/damage/isobrittle_generic.yaml index 98c9b9763..95c0e8b61 100644 --- a/examples/config/phase/damage/isobrittle_generic.yaml +++ b/examples/config/phase/damage/isobrittle_generic.yaml @@ -1,9 +1,7 @@ type: isobrittle W_crit: 1400000.0 -m: 1.0 -isoBrittle_atol: 0.01 output: [f_phi] -K_11: 1.0 +D_11: 1.0 mu: 0.001 diff --git a/python/damask/VERSION b/python/damask/VERSION index 862500995..026bbcf9f 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha4-94-g63fee141b +v3.0.0-alpha4-137-gb69b85754 diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 213b56a54..aa532859a 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -189,9 +189,11 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_eye(6) else validCalculation - if (debugCPFEM%extensive) & - print'(a,i8,1x,i2)', '<< CPFEM >> calculation for elFE ip ',elFE,ip - call materialpoint_stressAndItsTangent(dt,[ip,ip],[elCP,elCP]) + if (debugCPFEM%extensive) print'(a,i8,1x,i2)', '<< CPFEM >> calculation for elFE ip ',elFE,ip + call homogenization_mechanical_response(dt,[ip,ip],[elCP,elCP]) + if (.not. terminallyIll) & + call homogenization_mechanical_response2(dt,[ip,ip],[elCP,elCP]) + terminalIllness: if (terminallyIll) then diff --git a/src/IO.f90 b/src/IO.f90 index 3fbeb4300..399e2e6df 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -81,12 +81,7 @@ function IO_readlines(fileName) result(fileContent) rawData = IO_read(fileName) -!-------------------------------------------------------------------------------------------------- -! count lines to allocate string array - N_lines = 0 - do l=1, len(rawData) - if (rawData(l:l) == IO_EOL) N_lines = N_lines+1 - enddo + N_lines = count([(rawData(l:l) == IO_EOL,l=1,len(rawData))]) allocate(fileContent(N_lines)) !-------------------------------------------------------------------------------------------------- @@ -261,7 +256,7 @@ pure function IO_lc(string) do i=1,len(string) n = index(UPPER,string(i:i)) - if(n/=0) then + if (n/=0) then IO_lc(i:i) = LOWER(n:n) else IO_lc(i:i) = string(i:i) diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index e98b2d818..9c31b6f26 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -482,20 +482,13 @@ subroutine getMaskedTensor(values,mask,tensor) values = 0.0 - if (tensor%length == 9) then ! temporary support for deprecated 1D tensor - do i = 1,9 - mask((i-1)/3+1,mod(i-1,3)+1) = tensor%get_asString(i) /= 'x' - if (mask((i-1)/3+1,mod(i-1,3)+1)) values((i-1)/3+1,mod(i-1,3)+1) = tensor%get_asFloat(i) + do i = 1,3 + row => tensor%get(i) + do j = 1,3 + mask(i,j) = row%get_asString(j) /= 'x' + if (mask(i,j)) values(i,j) = row%get_asFloat(j) enddo - else - do i = 1,3 - row => tensor%get(i) - do j = 1,3 - mask(i,j) = row%get_asString(j) /= 'x' ! ToDo change to np.masked behavior - if (mask(i,j)) values(i,j) = row%get_asFloat(j) - enddo - enddo - endif + enddo end subroutine diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 77678137d..111e9fece 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -261,7 +261,7 @@ subroutine grid_mechanical_FEM_init F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) endif restartRead - homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent + homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for homogenization_mechanical_response call utilities_updateCoords(F) call utilities_constitutiveResponse(P_current,P_av,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2 F, & ! target F @@ -490,8 +490,8 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,i BCTol err_div = fnorm*sqrt(wgt)*geomSize(1)/scaledGeomSize(1)/detJ - divTol = max(maxval(abs(P_av))*num%eps_div_rtol ,num%eps_div_atol) - BCTol = max(maxval(abs(P_av))*num%eps_stress_rtol,num%eps_stress_atol) + divTol = max(maxval(abs(P_av))*num%eps_div_rtol, num%eps_div_atol) + BCTol = max(maxval(abs(P_av))*num%eps_stress_rtol, num%eps_stress_atol) if ((totalIter >= num%itmin .and. all([err_div/divTol, err_BC/BCTol] < 1.0_pReal)) & .or. terminallyIll) then @@ -502,8 +502,6 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,i reason = 0 endif -!-------------------------------------------------------------------------------------------------- -! report print'(1/,a)', ' ... reporting .............................................................' print'(1/,a,f12.2,a,es8.2,a,es9.2,a)', ' error divergence = ', & err_div/divTol, ' (',err_div,' / m, tol = ',divTol,')' diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 83b961023..c4fce61a7 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -211,7 +211,7 @@ subroutine grid_mechanical_spectral_basic_init F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) endif restartRead - homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent + homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for homogenization_mechanical_response call utilities_updateCoords(reshape(F,shape(F_lastInc))) call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 reshape(F,shape(F_lastInc)), & ! target F @@ -429,8 +429,8 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm divTol, & BCTol - divTol = max(maxval(abs(P_av))*num%eps_div_rtol ,num%eps_div_atol) - BCTol = max(maxval(abs(P_av))*num%eps_stress_rtol,num%eps_stress_atol) + divTol = max(maxval(abs(P_av))*num%eps_div_rtol, num%eps_div_atol) + BCTol = max(maxval(abs(P_av))*num%eps_stress_rtol, num%eps_stress_atol) if ((totalIter >= num%itmin .and. all([err_div/divTol, err_BC/BCTol] < 1.0_pReal)) & .or. terminallyIll) then @@ -441,8 +441,6 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm reason = 0 endif -!-------------------------------------------------------------------------------------------------- -! report print'(1/,a)', ' ... reporting .............................................................' print'(1/,a,f12.2,a,es8.2,a,es9.2,a)', ' error divergence = ', & err_div/divTol, ' (',err_div,' / m, tol = ',divTol,')' diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 07edc84b5..a52554e3c 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -237,7 +237,7 @@ subroutine grid_mechanical_spectral_polarisation_init F_tau_lastInc = 2.0_pReal*F_lastInc endif restartRead - homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent + homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for homogenization_mechanical_response call utilities_updateCoords(reshape(F,shape(F_lastInc))) call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 reshape(F,shape(F_lastInc)), & ! target F @@ -487,9 +487,9 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm divTol, & BCTol - curlTol = max(maxval(abs(F_aim-math_I3))*num%eps_curl_rtol ,num%eps_curl_atol) - divTol = max(maxval(abs(P_av)) *num%eps_div_rtol ,num%eps_div_atol) - BCTol = max(maxval(abs(P_av)) *num%eps_stress_rtol,num%eps_stress_atol) + curlTol = max(maxval(abs(F_aim-math_I3))*num%eps_curl_rtol, num%eps_curl_atol) + divTol = max(maxval(abs(P_av))*num%eps_div_rtol, num%eps_div_atol) + BCTol = max(maxval(abs(P_av))*num%eps_stress_rtol, num%eps_stress_atol) if ((totalIter >= num%itmin .and. all([err_div/divTol, err_curl/curlTol, err_BC/BCTol] < 1.0_pReal)) & .or. terminallyIll) then @@ -500,8 +500,6 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm reason = 0 endif -!-------------------------------------------------------------------------------------------------- -! report print'(1/,a)', ' ... reporting .............................................................' print'(1/,a,f12.2,a,es8.2,a,es9.2,a)', ' error divergence = ', & err_div/divTol, ' (',err_div, ' / m, tol = ',divTol,')' @@ -542,13 +540,13 @@ subroutine formResidual(in, FandF_tau, & !--------------------------------------------------------------------------------------------------- F => FandF_tau(1:3,1:3,1,& - XG_RANGE,YG_RANGE,ZG_RANGE) + XG_RANGE,YG_RANGE,ZG_RANGE) F_tau => FandF_tau(1:3,1:3,2,& - XG_RANGE,YG_RANGE,ZG_RANGE) + XG_RANGE,YG_RANGE,ZG_RANGE) residual_F => residuum(1:3,1:3,1,& - X_RANGE, Y_RANGE, Z_RANGE) + X_RANGE, Y_RANGE, Z_RANGE) residual_F_tau => residuum(1:3,1:3,2,& - X_RANGE, Y_RANGE, Z_RANGE) + X_RANGE, Y_RANGE, Z_RANGE) F_av = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr) diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 061cc9d74..bd37e2654 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -815,10 +815,14 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& homogenization_F = reshape(F,[3,3,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field - call materialpoint_stressAndItsTangent(timeinc,[1,1],[1,product(grid(1:2))*grid3]) ! calculate P field + call homogenization_mechanical_response(timeinc,[1,1],[1,product(grid(1:2))*grid3]) ! calculate P field + if (.not. terminallyIll) & + call homogenization_thermal_response(timeinc,[1,1],[1,product(grid(1:2))*grid3]) + if (.not. terminallyIll) & + call homogenization_mechanical_response2(timeinc,[1,1],[1,product(grid(1:2))*grid3]) P = reshape(homogenization_P, [3,3,grid(1),grid(2),grid3]) - P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt ! average of P + P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr) if (debugRotation) print'(/,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', & ' Piola--Kirchhoff stress (lab) / MPa =', transpose(P_av)*1.e-6_pReal diff --git a/src/homogenization.f90 b/src/homogenization.f90 index f34834272..924e45989 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -41,15 +41,6 @@ module homogenization integer(kind(DAMAGE_none_ID)), dimension(:), allocatable :: & damage_type !< nonlocal damage model - type, private :: tNumerics_damage - real(pReal) :: & - charLength !< characteristic length scale for gradient problems - end type tNumerics_damage - - type(tNumerics_damage), private :: & - num_damage - - logical, public :: & terminallyIll = .false. !< at least one material point is terminally ill @@ -101,8 +92,8 @@ module homogenization integer, intent(in) :: ce end subroutine damage_partition - module subroutine mechanical_homogenize(dt,ce) - real(pReal), intent(in) :: dt + module subroutine mechanical_homogenize(Delta_t,ce) + real(pReal), intent(in) :: Delta_t integer, intent(in) :: & ce !< cell end subroutine mechanical_homogenize @@ -178,7 +169,9 @@ module homogenization public :: & homogenization_init, & - materialpoint_stressAndItsTangent, & + homogenization_mechanical_response, & + homogenization_mechanical_response2, & + homogenization_thermal_response, & homogenization_mu_T, & homogenization_K_T, & homogenization_f_T, & @@ -227,24 +220,24 @@ end subroutine homogenization_init !-------------------------------------------------------------------------------------------------- -!> @brief parallelized calculation of stress and corresponding tangent at material points +!> @brief !-------------------------------------------------------------------------------------------------- -subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execElem) +subroutine homogenization_mechanical_response(Delta_t,FEsolving_execIP,FEsolving_execElem) - real(pReal), intent(in) :: dt !< time increment + real(pReal), intent(in) :: Delta_t !< time increment integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP integer :: & NiterationMPstate, & ip, & !< integration point number el, & !< element number - co, ce, ho, en, ph + co, ce, ho, en logical :: & converged logical, dimension(2) :: & doneAndHappy - !$OMP PARALLEL - !$OMP DO PRIVATE(ce,en,ho,NiterationMPstate,converged,doneAndHappy) + + !$OMP PARALLEL DO PRIVATE(ce,en,ho,NiterationMPstate,converged,doneAndHappy) do el = FEsolving_execElem(1),FEsolving_execElem(2) do ip = FEsolving_execIP(1),FEsolving_execIP(2) @@ -266,65 +259,89 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE NiterationMPstate = NiterationMPstate + 1 call mechanical_partition(homogenization_F(1:3,1:3,ce),ce) - converged = .true. - do co = 1, homogenization_Nconstituents(ho) - converged = converged .and. crystallite_stress(dt,co,ip,el) - enddo - + converged = all([(phase_mechanical_constitutive(Delta_t,co,ip,el),co=1,homogenization_Nconstituents(ho))]) if (converged) then - doneAndHappy = mechanical_updateState(dt,homogenization_F(1:3,1:3,ce),ce) + doneAndHappy = mechanical_updateState(Delta_t,homogenization_F(1:3,1:3,ce),ce) converged = all(doneAndHappy) else doneAndHappy = [.true.,.false.] endif - enddo convergenceLooping + + converged = converged .and. all([(phase_damage_constitutive(Delta_t,co,ip,el),co=1,homogenization_Nconstituents(ho))]) + if (.not. converged) then - if (.not. terminallyIll) print*, ' Integration point ', ip,' at element ', el, ' terminally ill' + if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill' terminallyIll = .true. endif enddo enddo - !$OMP END DO + !$OMP END PARALLEL DO - if (.not. terminallyIll) then - !$OMP DO PRIVATE(ho,ph,ce) - do el = FEsolving_execElem(1),FEsolving_execElem(2) - if (terminallyIll) continue - do ip = FEsolving_execIP(1),FEsolving_execIP(2) - ce = (el-1)*discretization_nIPs + ip - ho = material_homogenizationID(ce) - call thermal_partition(ce) - do co = 1, homogenization_Nconstituents(ho) - ph = material_phaseID(co,ce) - if (.not. thermal_stress(dt,ph,material_phaseMemberAt(co,ip,el))) then - if (.not. terminallyIll) & ! so first signals terminally ill... - print*, ' Integration point ', ip,' at element ', el, ' terminally ill' - terminallyIll = .true. ! ...and kills all others - endif - enddo +end subroutine homogenization_mechanical_response + + +!-------------------------------------------------------------------------------------------------- +!> @brief +!-------------------------------------------------------------------------------------------------- +subroutine homogenization_thermal_response(Delta_t,FEsolving_execIP,FEsolving_execElem) + + real(pReal), intent(in) :: Delta_t !< time increment + integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP + integer :: & + ip, & !< integration point number + el, & !< element number + co, ce, ho + + + !$OMP PARALLEL DO PRIVATE(ho,ce) + do el = FEsolving_execElem(1),FEsolving_execElem(2) + if (terminallyIll) continue + do ip = FEsolving_execIP(1),FEsolving_execIP(2) + ce = (el-1)*discretization_nIPs + ip + ho = material_homogenizationID(ce) + call thermal_partition(ce) + do co = 1, homogenization_Nconstituents(ho) + if (.not. phase_thermal_constitutive(Delta_t,material_phaseID(co,ce),material_phaseEntry(co,ce))) then + if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill' + terminallyIll = .true. + endif enddo enddo - !$OMP END DO + enddo + !$OMP END PARALLEL DO - !$OMP DO PRIVATE(ho,ce) - elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2) - IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2) - ce = (el-1)*discretization_nIPs + ip - ho = material_homogenizationID(ce) - do co = 1, homogenization_Nconstituents(ho) - call crystallite_orientations(co,ip,el) - enddo - call mechanical_homogenize(dt,ce) - enddo IpLooping3 - enddo elementLooping3 - !$OMP END DO - else - print'(/,a,/)', ' << HOMOG >> Material Point terminally ill' - endif - !$OMP END PARALLEL +end subroutine homogenization_thermal_response -end subroutine materialpoint_stressAndItsTangent + +!-------------------------------------------------------------------------------------------------- +!> @brief +!-------------------------------------------------------------------------------------------------- +subroutine homogenization_mechanical_response2(Delta_t,FEsolving_execIP,FEsolving_execElem) + + real(pReal), intent(in) :: Delta_t !< time increment + integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP + integer :: & + ip, & !< integration point number + el, & !< element number + co, ce, ho + + + !$OMP PARALLEL DO PRIVATE(ho,ce) + elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2) + IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2) + ce = (el-1)*discretization_nIPs + ip + ho = material_homogenizationID(ce) + do co = 1, homogenization_Nconstituents(ho) + call crystallite_orientations(co,ip,el) + enddo + call mechanical_homogenize(Delta_t,ce) + enddo IpLooping3 + enddo elementLooping3 + !$OMP END PARALLEL DO + + +end subroutine homogenization_mechanical_response2 !-------------------------------------------------------------------------------------------------- diff --git a/src/homogenization_damage.f90 b/src/homogenization_damage.f90 index f95322bfe..251cc72dc 100644 --- a/src/homogenization_damage.f90 +++ b/src/homogenization_damage.f90 @@ -37,8 +37,7 @@ module subroutine damage_init() class(tNode), pointer :: & configHomogenizations, & configHomogenization, & - configHomogenizationDamage, & - num_generic + configHomogenizationDamage integer :: ho,Nmembers @@ -70,11 +69,6 @@ module subroutine damage_init() end associate enddo -!------------------------------------------------------------------------------------ -! read numerics parameter - num_generic => config_numerics%get('generic',defaultVal= emptyDict) - num_damage%charLength = num_generic%get_asFloat('charLength',defaultVal=1.0_pReal) - call pass_init() end subroutine damage_init @@ -119,8 +113,7 @@ module function homogenization_K_phi(ce) result(K) real(pReal), dimension(3,3) :: K - K = phase_K_phi(1,ce) & - * num_damage%charLength**2 + K = phase_K_phi(1,ce) end function homogenization_K_phi diff --git a/src/homogenization_mechanical.f90 b/src/homogenization_mechanical.f90 index 7ea1f74e9..0138d2468 100644 --- a/src/homogenization_mechanical.f90 +++ b/src/homogenization_mechanical.f90 @@ -123,21 +123,21 @@ end subroutine mechanical_partition !-------------------------------------------------------------------------------------------------- !> @brief Average P and dPdF from the individual constituents. !-------------------------------------------------------------------------------------------------- -module subroutine mechanical_homogenize(dt,ce) +module subroutine mechanical_homogenize(Delta_t,ce) - real(pReal), intent(in) :: dt + real(pReal), intent(in) :: Delta_t integer, intent(in) :: ce integer :: co homogenization_P(1:3,1:3,ce) = phase_P(1,ce) - homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = phase_mechanical_dPdF(dt,1,ce) + homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = phase_mechanical_dPdF(Delta_t,1,ce) do co = 2, homogenization_Nconstituents(material_homogenizationID(ce)) homogenization_P(1:3,1:3,ce) = homogenization_P(1:3,1:3,ce) & + phase_P(co,ce) homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = homogenization_dPdF(1:3,1:3,1:3,1:3,ce) & - + phase_mechanical_dPdF(dt,co,ce) + + phase_mechanical_dPdF(Delta_t,co,ce) enddo homogenization_P(1:3,1:3,ce) = homogenization_P(1:3,1:3,ce) & diff --git a/src/lattice.f90 b/src/lattice.f90 index d366674c7..725c36ba8 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -3,7 +3,7 @@ !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH -!> @brief contains lattice structure definitions including Schmid matrices for slip, twin, trans, +!> @brief contains lattice definitions including Schmid matrices for slip, twin, trans, ! and cleavage as well as interaction among the various systems !-------------------------------------------------------------------------------------------------- module lattice @@ -365,12 +365,6 @@ module lattice 1, 1, 1, 1,-2, 1 & ],pReal),shape(BCT_SYSTEMSLIP)) !< bct slip systems for c/a = 0.5456 (Sn), sorted by Bieler 2009 (https://doi.org/10.1007/s11664-009-0909-x) -! SHOULD NOT BE PART OF LATTICE BEGIN - real(pReal), dimension(:), allocatable, public, protected :: & - lattice_mu, lattice_nu - real(pReal), dimension(:,:,:), allocatable, public, protected :: & - lattice_C66 -! SHOULD NOT BE PART OF LATTICE END interface lattice_forestProjection_edge module procedure slipProjection_transverse @@ -384,7 +378,8 @@ module lattice lattice_init, & lattice_equivalent_nu, & lattice_equivalent_mu, & - lattice_applyLatticeSymmetry33, & + lattice_symmetrize_33, & + lattice_symmetrize_C66, & lattice_SchmidMatrix_slip, & lattice_SchmidMatrix_twin, & lattice_SchmidMatrix_trans, & @@ -414,53 +409,9 @@ contains !-------------------------------------------------------------------------------------------------- subroutine lattice_init - integer :: Nphases, ph,i - class(tNode), pointer :: & - phases, & - phase, & - mech, & - elasticity - print'(/,a)', ' <<<+- lattice init -+>>>'; flush(IO_STDOUT) -! SHOULD NOT BE PART OF LATTICE BEGIN - - phases => config_material%get('phase') - Nphases = phases%length - - allocate(lattice_C66(6,6,Nphases), source=0.0_pReal) - - allocate(lattice_mu, lattice_nu,& - source=[(0.0_pReal,i=1,Nphases)]) - - do ph = 1, phases%length - phase => phases%get(ph) - mech => phase%get('mechanical') - elasticity => mech%get('elastic') - lattice_C66(1,1,ph) = elasticity%get_asFloat('C_11') - lattice_C66(1,2,ph) = elasticity%get_asFloat('C_12') - lattice_C66(4,4,ph) = elasticity%get_asFloat('C_44') - - lattice_C66(1,3,ph) = elasticity%get_asFloat('C_13',defaultVal=0.0_pReal) - lattice_C66(2,3,ph) = elasticity%get_asFloat('C_23',defaultVal=0.0_pReal) - lattice_C66(3,3,ph) = elasticity%get_asFloat('C_33',defaultVal=0.0_pReal) - lattice_C66(6,6,ph) = elasticity%get_asFloat('C_66',defaultVal=0.0_pReal) - - lattice_C66(1:6,1:6,ph) = applyLatticeSymmetryC66(lattice_C66(1:6,1:6,ph),phase%get_asString('lattice')) - - lattice_nu(ph) = lattice_equivalent_nu(lattice_C66(1:6,1:6,ph),'voigt') - lattice_mu(ph) = lattice_equivalent_mu(lattice_C66(1:6,1:6,ph),'voigt') - - lattice_C66(1:6,1:6,ph) = math_sym3333to66(math_Voigt66to3333(lattice_C66(1:6,1:6,ph))) ! Literature data is in Voigt notation - do i = 1, 6 - if (abs(lattice_C66(i,i,ph)) @brief Characteristic shear for twinning !-------------------------------------------------------------------------------------------------- -function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(characteristicShear) +function lattice_characteristicShear_Twin(Ntwin,lattice,CoverA) result(characteristicShear) integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family - character(len=*), intent(in) :: structure !< lattice structure + character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), dimension(sum(Ntwin)) :: characteristicShear @@ -513,7 +464,7 @@ function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(charact myFamilies: do f = 1,size(Ntwin,1) mySystems: do s = 1,Ntwin(f) a = a + 1 - select case(structure) + select case(lattice) case('cF','cI') characteristicShear(a) = 0.5_pReal*sqrt(2.0_pReal) case('hP') @@ -531,7 +482,7 @@ function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(charact characteristicShear(a) = 2.0_pReal*(cOverA**2.0_pReal-2.0_pReal)/3.0_pReal/cOverA end select case default - call IO_error(137,ext_msg='lattice_characteristicShear_Twin: '//trim(structure)) + call IO_error(137,ext_msg='lattice_characteristicShear_Twin: '//trim(lattice)) end select enddo mySystems enddo myFamilies @@ -542,10 +493,10 @@ end function lattice_characteristicShear_Twin !-------------------------------------------------------------------------------------------------- !> @brief Rotated elasticity matrices for twinning in 66-vector notation !-------------------------------------------------------------------------------------------------- -function lattice_C66_twin(Ntwin,C66,structure,CoverA) +function lattice_C66_twin(Ntwin,C66,lattice,CoverA) integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family - character(len=*), intent(in) :: structure !< lattice structure + character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) real(pReal), dimension(6,6), intent(in) :: C66 !< unrotated parent stiffness matrix real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), dimension(6,6,sum(Ntwin)) :: lattice_C66_twin @@ -554,18 +505,18 @@ function lattice_C66_twin(Ntwin,C66,structure,CoverA) type(rotation) :: R integer :: i - select case(structure) + select case(lattice) case('cF') coordinateSystem = buildCoordinateSystem(Ntwin,FCC_NSLIPSYSTEM,FCC_SYSTEMTWIN,& - structure,0.0_pReal) + lattice,0.0_pReal) case('cI') coordinateSystem = buildCoordinateSystem(Ntwin,BCC_NSLIPSYSTEM,BCC_SYSTEMTWIN,& - structure,0.0_pReal) + lattice,0.0_pReal) case('hP') coordinateSystem = buildCoordinateSystem(Ntwin,HEX_NSLIPSYSTEM,HEX_SYSTEMTWIN,& - structure,cOverA) + lattice,cOverA) case default - call IO_error(137,ext_msg='lattice_C66_twin: '//trim(structure)) + call IO_error(137,ext_msg='lattice_C66_twin: '//trim(lattice)) end select do i = 1, sum(Ntwin) @@ -579,11 +530,11 @@ end function lattice_C66_twin !-------------------------------------------------------------------------------------------------- !> @brief Rotated elasticity matrices for transformation in 66-vector notation !-------------------------------------------------------------------------------------------------- -function lattice_C66_trans(Ntrans,C_parent66,structure_target, & +function lattice_C66_trans(Ntrans,C_parent66,lattice_target, & cOverA_trans,a_bcc,a_fcc) integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family - character(len=*), intent(in) :: structure_target !< lattice structure + character(len=2), intent(in) :: lattice_target !< Bravais lattice (Pearson symbol) real(pReal), dimension(6,6), intent(in) :: C_parent66 real(pReal), dimension(6,6,sum(Ntrans)) :: lattice_C66_trans @@ -595,9 +546,9 @@ function lattice_C66_trans(Ntrans,C_parent66,structure_target, & !-------------------------------------------------------------------------------------------------- ! elasticity matrix of the target phase in cube orientation - if (structure_target == 'hP') then + if (lattice_target == 'hP') then if (cOverA_trans < 1.0_pReal .or. cOverA_trans > 2.0_pReal) & - call IO_error(131,ext_msg='lattice_C66_trans: '//trim(structure_target)) + call IO_error(131,ext_msg='lattice_C66_trans: '//trim(lattice_target)) C_bar66(1,1) = (C_parent66(1,1) + C_parent66(1,2) + 2.0_pReal*C_parent66(4,4))/2.0_pReal C_bar66(1,2) = (C_parent66(1,1) + 5.0_pReal*C_parent66(1,2) - 2.0_pReal*C_parent66(4,4))/6.0_pReal C_bar66(3,3) = (C_parent66(1,1) + 2.0_pReal*C_parent66(1,2) + 4.0_pReal*C_parent66(4,4))/3.0_pReal @@ -611,13 +562,13 @@ function lattice_C66_trans(Ntrans,C_parent66,structure_target, & C_target_unrotated66(1,3) = C_bar66(1,3) C_target_unrotated66(3,3) = C_bar66(3,3) C_target_unrotated66(4,4) = C_bar66(4,4) - C_bar66(1,4)**2.0_pReal/(0.5_pReal*(C_bar66(1,1) - C_bar66(1,2))) - C_target_unrotated66 = applyLatticeSymmetryC66(C_target_unrotated66,'hP') - elseif (structure_target == 'cI') then + C_target_unrotated66 = lattice_symmetrize_C66(C_target_unrotated66,'hP') + elseif (lattice_target == 'cI') then if (a_bcc <= 0.0_pReal .or. a_fcc <= 0.0_pReal) & - call IO_error(134,ext_msg='lattice_C66_trans: '//trim(structure_target)) + call IO_error(134,ext_msg='lattice_C66_trans: '//trim(lattice_target)) C_target_unrotated66 = C_parent66 else - call IO_error(137,ext_msg='lattice_C66_trans : '//trim(structure_target)) + call IO_error(137,ext_msg='lattice_C66_trans : '//trim(lattice_target)) endif do i = 1, 6 @@ -686,11 +637,11 @@ end function lattice_nonSchmidMatrix !> @brief Slip-slip interaction matrix !> details only active slip systems are considered !-------------------------------------------------------------------------------------------------- -function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) result(interactionMatrix) +function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(interactionMatrix) integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family real(pReal), dimension(:), intent(in) :: interactionValues !< values for slip-slip interaction - character(len=*), intent(in) :: structure !< lattice structure + character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) real(pReal), dimension(sum(Nslip),sum(Nslip)) :: interactionMatrix integer, dimension(:), allocatable :: NslipMax @@ -908,7 +859,7 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) resul ],shape(BCT_INTERACTIONSLIPSLIP)) - select case(structure) + select case(lattice) case('cF') interactionTypes = FCC_INTERACTIONSLIPSLIP NslipMax = FCC_NSLIPSYSTEM @@ -922,7 +873,7 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) resul interactionTypes = BCT_INTERACTIONSLIPSLIP NslipMax = BCT_NSLIPSYSTEM case default - call IO_error(137,ext_msg='lattice_interaction_SlipBySlip: '//trim(structure)) + call IO_error(137,ext_msg='lattice_interaction_SlipBySlip: '//trim(lattice)) end select interactionMatrix = buildInteraction(Nslip,Nslip,NslipMax,NslipMax,interactionValues,interactionTypes) @@ -934,11 +885,11 @@ end function lattice_interaction_SlipBySlip !> @brief Twin-twin interaction matrix !> details only active twin systems are considered !-------------------------------------------------------------------------------------------------- -function lattice_interaction_TwinByTwin(Ntwin,interactionValues,structure) result(interactionMatrix) +function lattice_interaction_TwinByTwin(Ntwin,interactionValues,lattice) result(interactionMatrix) integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family real(pReal), dimension(:), intent(in) :: interactionValues !< values for twin-twin interaction - character(len=*), intent(in) :: structure !< lattice structure + character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) real(pReal), dimension(sum(Ntwin),sum(Ntwin)) :: interactionMatrix integer, dimension(:), allocatable :: NtwinMax @@ -1009,7 +960,7 @@ function lattice_interaction_TwinByTwin(Ntwin,interactionValues,structure) resul 20,20,20,20,20,20, 19,19,19,19,19,19, 18,18,18,18,18,18, 17,17,17,17,17,16 & ],shape(HEX_INTERACTIONTWINTWIN)) !< Twin-twin interaction types for hex - select case(structure) + select case(lattice) case('cF') interactionTypes = FCC_INTERACTIONTWINTWIN NtwinMax = FCC_NTWINSYSTEM @@ -1020,7 +971,7 @@ function lattice_interaction_TwinByTwin(Ntwin,interactionValues,structure) resul interactionTypes = HEX_INTERACTIONTWINTWIN NtwinMax = HEX_NTWINSYSTEM case default - call IO_error(137,ext_msg='lattice_interaction_TwinByTwin: '//trim(structure)) + call IO_error(137,ext_msg='lattice_interaction_TwinByTwin: '//trim(lattice)) end select interactionMatrix = buildInteraction(Ntwin,Ntwin,NtwinMax,NtwinMax,interactionValues,interactionTypes) @@ -1032,11 +983,11 @@ end function lattice_interaction_TwinByTwin !> @brief Trans-trans interaction matrix !> details only active trans systems are considered !-------------------------------------------------------------------------------------------------- -function lattice_interaction_TransByTrans(Ntrans,interactionValues,structure) result(interactionMatrix) +function lattice_interaction_TransByTrans(Ntrans,interactionValues,lattice) result(interactionMatrix) integer, dimension(:), intent(in) :: Ntrans !< number of active trans systems per family real(pReal), dimension(:), intent(in) :: interactionValues !< values for trans-trans interaction - character(len=*), intent(in) :: structure !< lattice structure (parent crystal) + character(len=2), intent(in) :: lattice ! @brief Slip-twin interaction matrix !> details only active slip and twin systems are considered !-------------------------------------------------------------------------------------------------- -function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,structure) result(interactionMatrix) +function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,lattice) result(interactionMatrix) integer, dimension(:), intent(in) :: Nslip, & !< number of active slip systems per family Ntwin !< number of active twin systems per family real(pReal), dimension(:), intent(in) :: interactionValues !< values for slip-twin interaction - character(len=*), intent(in) :: structure !< lattice structure + character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) real(pReal), dimension(sum(Nslip),sum(Ntwin)) :: interactionMatrix integer, dimension(:), allocatable :: NslipMax, & @@ -1211,7 +1162,7 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,structure) 21,21,21,21,21,21, 22,22,22,22,22,22, 23,23,23,23,23,23, 24,24,24,24,24,24 & ],shape(HEX_INTERACTIONSLIPTWIN)) !< Slip-twin interaction types for hex - select case(structure) + select case(lattice) case('cF') interactionTypes = FCC_INTERACTIONSLIPTWIN NslipMax = FCC_NSLIPSYSTEM @@ -1225,7 +1176,7 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,structure) NslipMax = HEX_NSLIPSYSTEM NtwinMax = HEX_NTWINSYSTEM case default - call IO_error(137,ext_msg='lattice_interaction_SlipByTwin: '//trim(structure)) + call IO_error(137,ext_msg='lattice_interaction_SlipByTwin: '//trim(lattice)) end select interactionMatrix = buildInteraction(Nslip,Ntwin,NslipMax,NtwinMax,interactionValues,interactionTypes) @@ -1237,12 +1188,12 @@ end function lattice_interaction_SlipByTwin !> @brief Slip-trans interaction matrix !> details only active slip and trans systems are considered !-------------------------------------------------------------------------------------------------- -function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,structure) result(interactionMatrix) +function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,lattice) result(interactionMatrix) integer, dimension(:), intent(in) :: Nslip, & !< number of active slip systems per family Ntrans !< number of active trans systems per family real(pReal), dimension(:), intent(in) :: interactionValues !< values for slip-trans interaction - character(len=*), intent(in) :: structure !< lattice structure (parent crystal) + character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) (parent crystal) real(pReal), dimension(sum(Nslip),sum(Ntrans)) :: interactionMatrix integer, dimension(:), allocatable :: NslipMax, & @@ -1272,13 +1223,13 @@ function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,structur 4,4,4,4,4,4,4,4,4,4,4,4 & ],shape(FCC_INTERACTIONSLIPTRANS)) !< Slip-trans interaction types for fcc - select case(structure) + select case(lattice) case('cF') interactionTypes = FCC_INTERACTIONSLIPTRANS NslipMax = FCC_NSLIPSYSTEM NtransMax = FCC_NTRANSSYSTEM case default - call IO_error(137,ext_msg='lattice_interaction_SlipByTrans: '//trim(structure)) + call IO_error(137,ext_msg='lattice_interaction_SlipByTrans: '//trim(lattice)) end select interactionMatrix = buildInteraction(Nslip,Ntrans,NslipMax,NtransMax,interactionValues,interactionTypes) @@ -1290,12 +1241,12 @@ function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,structur !> @brief Twin-slip interaction matrix !> details only active twin and slip systems are considered !-------------------------------------------------------------------------------------------------- -function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,structure) result(interactionMatrix) +function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,lattice) result(interactionMatrix) integer, dimension(:), intent(in) :: Ntwin, & !< number of active twin systems per family Nslip !< number of active slip systems per family real(pReal), dimension(:), intent(in) :: interactionValues !< values for twin-twin interaction - character(len=*), intent(in) :: structure !< lattice structure + character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) real(pReal), dimension(sum(Ntwin),sum(Nslip)) :: interactionMatrix integer, dimension(:), allocatable :: NtwinMax, & @@ -1339,7 +1290,7 @@ function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,structure) 4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24 & ],shape(HEX_INTERACTIONTWINSLIP)) !< Twin-slip interaction types for hex - select case(structure) + select case(lattice) case('cF') interactionTypes = FCC_INTERACTIONTWINSLIP NtwinMax = FCC_NTWINSYSTEM @@ -1353,7 +1304,7 @@ function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,structure) NtwinMax = HEX_NTWINSYSTEM NslipMax = HEX_NSLIPSYSTEM case default - call IO_error(137,ext_msg='lattice_interaction_TwinBySlip: '//trim(structure)) + call IO_error(137,ext_msg='lattice_interaction_TwinBySlip: '//trim(lattice)) end select interactionMatrix = buildInteraction(Ntwin,Nslip,NtwinMax,NslipMax,interactionValues,interactionTypes) @@ -1365,10 +1316,10 @@ end function lattice_interaction_TwinBySlip !> @brief Schmid matrix for slip !> details only active slip systems are considered !-------------------------------------------------------------------------------------------------- -function lattice_SchmidMatrix_slip(Nslip,structure,cOverA) result(SchmidMatrix) +function lattice_SchmidMatrix_slip(Nslip,lattice,cOverA) result(SchmidMatrix) integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family - character(len=*), intent(in) :: structure !< lattice structure + character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) real(pReal), intent(in) :: cOverA real(pReal), dimension(3,3,sum(Nslip)) :: SchmidMatrix @@ -1377,7 +1328,7 @@ function lattice_SchmidMatrix_slip(Nslip,structure,cOverA) result(SchmidMatrix) integer, dimension(:), allocatable :: NslipMax integer :: i - select case(structure) + select case(lattice) case('cF') NslipMax = FCC_NSLIPSYSTEM slipSystems = FCC_SYSTEMSLIP @@ -1392,15 +1343,15 @@ function lattice_SchmidMatrix_slip(Nslip,structure,cOverA) result(SchmidMatrix) slipSystems = BCT_SYSTEMSLIP case default allocate(NslipMax(0)) - call IO_error(137,ext_msg='lattice_SchmidMatrix_slip: '//trim(structure)) + call IO_error(137,ext_msg='lattice_SchmidMatrix_slip: '//trim(lattice)) end select if (any(NslipMax(1:size(Nslip)) - Nslip < 0)) & - call IO_error(145,ext_msg='Nslip '//trim(structure)) + call IO_error(145,ext_msg='Nslip '//trim(lattice)) if (any(Nslip < 0)) & - call IO_error(144,ext_msg='Nslip '//trim(structure)) + call IO_error(144,ext_msg='Nslip '//trim(lattice)) - coordinateSystem = buildCoordinateSystem(Nslip,NslipMax,slipSystems,structure,cOverA) + coordinateSystem = buildCoordinateSystem(Nslip,NslipMax,slipSystems,lattice,cOverA) do i = 1, sum(Nslip) SchmidMatrix(1:3,1:3,i) = math_outer(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i)) @@ -1415,10 +1366,10 @@ end function lattice_SchmidMatrix_slip !> @brief Schmid matrix for twinning !> details only active twin systems are considered !-------------------------------------------------------------------------------------------------- -function lattice_SchmidMatrix_twin(Ntwin,structure,cOverA) result(SchmidMatrix) +function lattice_SchmidMatrix_twin(Ntwin,lattice,cOverA) result(SchmidMatrix) integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family - character(len=*), intent(in) :: structure !< lattice structure + character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), dimension(3,3,sum(Ntwin)) :: SchmidMatrix @@ -1427,7 +1378,7 @@ function lattice_SchmidMatrix_twin(Ntwin,structure,cOverA) result(SchmidMatrix) integer, dimension(:), allocatable :: NtwinMax integer :: i - select case(structure) + select case(lattice) case('cF') NtwinMax = FCC_NTWINSYSTEM twinSystems = FCC_SYSTEMTWIN @@ -1439,15 +1390,15 @@ function lattice_SchmidMatrix_twin(Ntwin,structure,cOverA) result(SchmidMatrix) twinSystems = HEX_SYSTEMTWIN case default allocate(NtwinMax(0)) - call IO_error(137,ext_msg='lattice_SchmidMatrix_twin: '//trim(structure)) + call IO_error(137,ext_msg='lattice_SchmidMatrix_twin: '//trim(lattice)) end select if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0)) & - call IO_error(145,ext_msg='Ntwin '//trim(structure)) + call IO_error(145,ext_msg='Ntwin '//trim(lattice)) if (any(Ntwin < 0)) & - call IO_error(144,ext_msg='Ntwin '//trim(structure)) + call IO_error(144,ext_msg='Ntwin '//trim(lattice)) - coordinateSystem = buildCoordinateSystem(Ntwin,NtwinMax,twinSystems,structure,cOverA) + coordinateSystem = buildCoordinateSystem(Ntwin,NtwinMax,twinSystems,lattice,cOverA) do i = 1, sum(Ntwin) SchmidMatrix(1:3,1:3,i) = math_outer(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i)) @@ -1462,24 +1413,24 @@ end function lattice_SchmidMatrix_twin !> @brief Schmid matrix for twinning !> details only active twin systems are considered !-------------------------------------------------------------------------------------------------- -function lattice_SchmidMatrix_trans(Ntrans,structure_target,cOverA,a_bcc,a_fcc) result(SchmidMatrix) +function lattice_SchmidMatrix_trans(Ntrans,lattice_target,cOverA,a_bcc,a_fcc) result(SchmidMatrix) integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family - character(len=*), intent(in) :: structure_target !< lattice structure + character(len=2), intent(in) :: lattice_target !< Bravais lattice (Pearson symbol) real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), dimension(3,3,sum(Ntrans)) :: SchmidMatrix real(pReal), dimension(3,3,sum(Ntrans)) :: devNull real(pReal) :: a_bcc, a_fcc - if (structure_target /= 'cI' .and. structure_target /= 'hP') & - call IO_error(137,ext_msg='lattice_SchmidMatrix_trans: '//trim(structure_target)) + if (lattice_target /= 'cI' .and. lattice_target /= 'hP') & + call IO_error(137,ext_msg='lattice_SchmidMatrix_trans: '//trim(lattice_target)) - if (structure_target == 'hP' .and. (cOverA < 1.0_pReal .or. cOverA > 2.0_pReal)) & - call IO_error(131,ext_msg='lattice_SchmidMatrix_trans: '//trim(structure_target)) + if (lattice_target == 'hP' .and. (cOverA < 1.0_pReal .or. cOverA > 2.0_pReal)) & + call IO_error(131,ext_msg='lattice_SchmidMatrix_trans: '//trim(lattice_target)) - if (structure_target == 'cI' .and. (a_bcc <= 0.0_pReal .or. a_fcc <= 0.0_pReal)) & - call IO_error(134,ext_msg='lattice_SchmidMatrix_trans: '//trim(structure_target)) + if (lattice_target == 'cI' .and. (a_bcc <= 0.0_pReal .or. a_fcc <= 0.0_pReal)) & + call IO_error(134,ext_msg='lattice_SchmidMatrix_trans: '//trim(lattice_target)) call buildTransformationSystem(devNull,SchmidMatrix,Ntrans,cOverA,a_fcc,a_bcc) @@ -1490,10 +1441,10 @@ end function lattice_SchmidMatrix_trans !> @brief Schmid matrix for cleavage !> details only active cleavage systems are considered !-------------------------------------------------------------------------------------------------- -function lattice_SchmidMatrix_cleavage(Ncleavage,structure,cOverA) result(SchmidMatrix) +function lattice_SchmidMatrix_cleavage(Ncleavage,lattice,cOverA) result(SchmidMatrix) integer, dimension(:), intent(in) :: Ncleavage !< number of active cleavage systems per family - character(len=*), intent(in) :: structure !< lattice structure + character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), dimension(3,3,3,sum(Ncleavage)) :: SchmidMatrix @@ -1502,7 +1453,7 @@ function lattice_SchmidMatrix_cleavage(Ncleavage,structure,cOverA) result(Schmid integer, dimension(:), allocatable :: NcleavageMax integer :: i - select case(structure) + select case(lattice) case('cF') NcleavageMax = FCC_NCLEAVAGESYSTEM cleavageSystems = FCC_SYSTEMCLEAVAGE @@ -1511,15 +1462,15 @@ function lattice_SchmidMatrix_cleavage(Ncleavage,structure,cOverA) result(Schmid cleavageSystems = BCC_SYSTEMCLEAVAGE case default allocate(NcleavageMax(0)) - call IO_error(137,ext_msg='lattice_SchmidMatrix_cleavage: '//trim(structure)) + call IO_error(137,ext_msg='lattice_SchmidMatrix_cleavage: '//trim(lattice)) end select if (any(NcleavageMax(1:size(Ncleavage)) - Ncleavage < 0)) & - call IO_error(145,ext_msg='Ncleavage '//trim(structure)) + call IO_error(145,ext_msg='Ncleavage '//trim(lattice)) if (any(Ncleavage < 0)) & - call IO_error(144,ext_msg='Ncleavage '//trim(structure)) + call IO_error(144,ext_msg='Ncleavage '//trim(lattice)) - coordinateSystem = buildCoordinateSystem(Ncleavage,NcleavageMax,cleavageSystems,structure,cOverA) + coordinateSystem = buildCoordinateSystem(Ncleavage,NcleavageMax,cleavageSystems,lattice,cOverA) do i = 1, sum(Ncleavage) SchmidMatrix(1:3,1:3,1,i) = math_outer(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i)) @@ -1533,16 +1484,16 @@ end function lattice_SchmidMatrix_cleavage !-------------------------------------------------------------------------------------------------- !> @brief Slip direction of slip systems (|| b) !-------------------------------------------------------------------------------------------------- -function lattice_slip_direction(Nslip,structure,cOverA) result(d) +function lattice_slip_direction(Nslip,lattice,cOverA) result(d) integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family - character(len=*), intent(in) :: structure !< lattice structure + character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), dimension(3,sum(Nslip)) :: d real(pReal), dimension(3,3,sum(Nslip)) :: coordinateSystem - coordinateSystem = coordinateSystem_slip(Nslip,structure,cOverA) + coordinateSystem = coordinateSystem_slip(Nslip,lattice,cOverA) d = coordinateSystem(1:3,1,1:sum(Nslip)) end function lattice_slip_direction @@ -1551,16 +1502,16 @@ end function lattice_slip_direction !-------------------------------------------------------------------------------------------------- !> @brief Normal direction of slip systems (|| n) !-------------------------------------------------------------------------------------------------- -function lattice_slip_normal(Nslip,structure,cOverA) result(n) +function lattice_slip_normal(Nslip,lattice,cOverA) result(n) integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family - character(len=*), intent(in) :: structure !< lattice structure + character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), dimension(3,sum(Nslip)) :: n real(pReal), dimension(3,3,sum(Nslip)) :: coordinateSystem - coordinateSystem = coordinateSystem_slip(Nslip,structure,cOverA) + coordinateSystem = coordinateSystem_slip(Nslip,lattice,cOverA) n = coordinateSystem(1:3,2,1:sum(Nslip)) end function lattice_slip_normal @@ -1569,16 +1520,16 @@ end function lattice_slip_normal !-------------------------------------------------------------------------------------------------- !> @brief Transverse direction of slip systems (|| t = b x n) !-------------------------------------------------------------------------------------------------- -function lattice_slip_transverse(Nslip,structure,cOverA) result(t) +function lattice_slip_transverse(Nslip,lattice,cOverA) result(t) integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family - character(len=*), intent(in) :: structure !< lattice structure + character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), dimension(3,sum(Nslip)) :: t real(pReal), dimension(3,3,sum(Nslip)) :: coordinateSystem - coordinateSystem = coordinateSystem_slip(Nslip,structure,cOverA) + coordinateSystem = coordinateSystem_slip(Nslip,lattice,cOverA) t = coordinateSystem(1:3,3,1:sum(Nslip)) end function lattice_slip_transverse @@ -1588,17 +1539,17 @@ end function lattice_slip_transverse !> @brief Labels for slip systems !> details only active slip systems are considered !-------------------------------------------------------------------------------------------------- -function lattice_labels_slip(Nslip,structure) result(labels) +function lattice_labels_slip(Nslip,lattice) result(labels) integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family - character(len=*), intent(in) :: structure !< lattice structure + character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=:), dimension(:), allocatable :: labels real(pReal), dimension(:,:), allocatable :: slipSystems integer, dimension(:), allocatable :: NslipMax - select case(structure) + select case(lattice) case('cF') NslipMax = FCC_NSLIPSYSTEM slipSystems = FCC_SYSTEMSLIP @@ -1612,13 +1563,13 @@ function lattice_labels_slip(Nslip,structure) result(labels) NslipMax = BCT_NSLIPSYSTEM slipSystems = BCT_SYSTEMSLIP case default - call IO_error(137,ext_msg='lattice_labels_slip: '//trim(structure)) + call IO_error(137,ext_msg='lattice_labels_slip: '//trim(lattice)) end select if (any(NslipMax(1:size(Nslip)) - Nslip < 0)) & - call IO_error(145,ext_msg='Nslip '//trim(structure)) + call IO_error(145,ext_msg='Nslip '//trim(lattice)) if (any(Nslip < 0)) & - call IO_error(144,ext_msg='Nslip '//trim(structure)) + call IO_error(144,ext_msg='Nslip '//trim(lattice)) labels = getLabels(Nslip,NslipMax,slipSystems) @@ -1626,19 +1577,19 @@ end function lattice_labels_slip !-------------------------------------------------------------------------------------------------- -!> @brief Return 3x3 tensor with symmetry according to given crystal structure +!> @brief Return 3x3 tensor with symmetry according to given Bravais lattice !-------------------------------------------------------------------------------------------------- -pure function lattice_applyLatticeSymmetry33(T,structure) result(T_sym) +pure function lattice_symmetrize_33(T,lattice) result(T_sym) real(pReal), dimension(3,3) :: T_sym real(pReal), dimension(3,3), intent(in) :: T - character(len=*), intent(in) :: structure + character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) T_sym = 0.0_pReal - select case(structure) + select case(lattice) case('cF','cI') T_sym(1,1) = T(1,1) T_sym(2,2) = T(1,1) @@ -1649,26 +1600,26 @@ pure function lattice_applyLatticeSymmetry33(T,structure) result(T_sym) T_sym(3,3) = T(3,3) end select -end function lattice_applyLatticeSymmetry33 +end function lattice_symmetrize_33 !-------------------------------------------------------------------------------------------------- -!> @brief Return stiffness matrix in 6x6 notation with symmetry according to given crystal structure +!> @brief Return stiffness matrix in 6x6 notation with symmetry according to given Bravais lattice !> @details J. A. Rayne and B. S. Chandrasekhar Phys. Rev. 120, 1658 Erratum Phys. Rev. 122, 1962 !-------------------------------------------------------------------------------------------------- -pure function applyLatticeSymmetryC66(C66,structure) result(C66_sym) +pure function lattice_symmetrize_C66(C66,lattice) result(C66_sym) real(pReal), dimension(6,6) :: C66_sym real(pReal), dimension(6,6), intent(in) :: C66 - character(len=*), intent(in) :: structure + character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) integer :: i,j C66_sym = 0.0_pReal - select case(structure) + select case(lattice) case ('cF','cI') C66_sym(1,1) = C66(1,1); C66_sym(2,2) = C66(1,1); C66_sym(3,3) = C66(1,1) C66_sym(1,2) = C66(1,2); C66_sym(1,3) = C66(1,2); C66_sym(2,3) = C66(1,2) @@ -1695,24 +1646,24 @@ pure function applyLatticeSymmetryC66(C66,structure) result(C66_sym) enddo enddo -end function applyLatticeSymmetryC66 +end function lattice_symmetrize_C66 !-------------------------------------------------------------------------------------------------- !> @brief Labels for twin systems !> details only active twin systems are considered !-------------------------------------------------------------------------------------------------- -function lattice_labels_twin(Ntwin,structure) result(labels) +function lattice_labels_twin(Ntwin,lattice) result(labels) integer, dimension(:), intent(in) :: Ntwin !< number of active slip systems per family - character(len=*), intent(in) :: structure !< lattice structure + character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=:), dimension(:), allocatable :: labels real(pReal), dimension(:,:), allocatable :: twinSystems integer, dimension(:), allocatable :: NtwinMax - select case(structure) + select case(lattice) case('cF') NtwinMax = FCC_NTWINSYSTEM twinSystems = FCC_SYSTEMTWIN @@ -1723,13 +1674,13 @@ function lattice_labels_twin(Ntwin,structure) result(labels) NtwinMax = HEX_NTWINSYSTEM twinSystems = HEX_SYSTEMTWIN case default - call IO_error(137,ext_msg='lattice_labels_twin: '//trim(structure)) + call IO_error(137,ext_msg='lattice_labels_twin: '//trim(lattice)) end select if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0)) & - call IO_error(145,ext_msg='Ntwin '//trim(structure)) + call IO_error(145,ext_msg='Ntwin '//trim(lattice)) if (any(Ntwin < 0)) & - call IO_error(144,ext_msg='Ntwin '//trim(structure)) + call IO_error(144,ext_msg='Ntwin '//trim(lattice)) labels = getLabels(Ntwin,NtwinMax,twinSystems) @@ -1740,18 +1691,18 @@ end function lattice_labels_twin !> @brief Projection of the transverse direction onto the slip plane !> @details: This projection is used to calculate forest hardening for edge dislocations !-------------------------------------------------------------------------------------------------- -function slipProjection_transverse(Nslip,structure,cOverA) result(projection) +function slipProjection_transverse(Nslip,lattice,cOverA) result(projection) integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family - character(len=*), intent(in) :: structure !< lattice structure + character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), dimension(sum(Nslip),sum(Nslip)) :: projection real(pReal), dimension(3,sum(Nslip)) :: n, t integer :: i, j - n = lattice_slip_normal (Nslip,structure,cOverA) - t = lattice_slip_transverse(Nslip,structure,cOverA) + n = lattice_slip_normal (Nslip,lattice,cOverA) + t = lattice_slip_transverse(Nslip,lattice,cOverA) do i=1, sum(Nslip); do j=1, sum(Nslip) projection(i,j) = abs(math_inner(n(:,i),t(:,j))) @@ -1764,18 +1715,18 @@ end function slipProjection_transverse !> @brief Projection of the slip direction onto the slip plane !> @details: This projection is used to calculate forest hardening for screw dislocations !-------------------------------------------------------------------------------------------------- -function slipProjection_direction(Nslip,structure,cOverA) result(projection) +function slipProjection_direction(Nslip,lattice,cOverA) result(projection) integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family - character(len=*), intent(in) :: structure !< lattice structure + character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), dimension(sum(Nslip),sum(Nslip)) :: projection real(pReal), dimension(3,sum(Nslip)) :: n, d integer :: i, j - n = lattice_slip_normal (Nslip,structure,cOverA) - d = lattice_slip_direction(Nslip,structure,cOverA) + n = lattice_slip_normal (Nslip,lattice,cOverA) + d = lattice_slip_direction(Nslip,lattice,cOverA) do i=1, sum(Nslip); do j=1, sum(Nslip) projection(i,j) = abs(math_inner(n(:,i),d(:,j))) @@ -1788,17 +1739,17 @@ end function slipProjection_direction !> @brief build a local coordinate system on slip systems !> @details Order: Direction, plane (normal), and common perpendicular !-------------------------------------------------------------------------------------------------- -function coordinateSystem_slip(Nslip,structure,cOverA) result(coordinateSystem) +function coordinateSystem_slip(Nslip,lattice,cOverA) result(coordinateSystem) integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family - character(len=*), intent(in) :: structure !< lattice structure + character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), dimension(3,3,sum(Nslip)) :: coordinateSystem real(pReal), dimension(:,:), allocatable :: slipSystems integer, dimension(:), allocatable :: NslipMax - select case(structure) + select case(lattice) case('cF') NslipMax = FCC_NSLIPSYSTEM slipSystems = FCC_SYSTEMSLIP @@ -1813,15 +1764,15 @@ function coordinateSystem_slip(Nslip,structure,cOverA) result(coordinateSystem) slipSystems = BCT_SYSTEMSLIP case default allocate(NslipMax(0)) - call IO_error(137,ext_msg='coordinateSystem_slip: '//trim(structure)) + call IO_error(137,ext_msg='coordinateSystem_slip: '//trim(lattice)) end select if (any(NslipMax(1:size(Nslip)) - Nslip < 0)) & - call IO_error(145,ext_msg='Nslip '//trim(structure)) + call IO_error(145,ext_msg='Nslip '//trim(lattice)) if (any(Nslip < 0)) & - call IO_error(144,ext_msg='Nslip '//trim(structure)) + call IO_error(144,ext_msg='Nslip '//trim(lattice)) - coordinateSystem = buildCoordinateSystem(Nslip,NslipMax,slipSystems,structure,cOverA) + coordinateSystem = buildCoordinateSystem(Nslip,NslipMax,slipSystems,lattice,cOverA) end function coordinateSystem_slip @@ -1873,15 +1824,15 @@ end function buildInteraction !> @brief Build a local coordinate system on slip, twin, trans, cleavage systems !> @details Order: Direction, plane (normal), and common perpendicular !-------------------------------------------------------------------------------------------------- -function buildCoordinateSystem(active,potential,system,structure,cOverA) +function buildCoordinateSystem(active,potential,system,lattice,cOverA) integer, dimension(:), intent(in) :: & active, & !< # of active systems per family potential !< # of potential systems per family real(pReal), dimension(:,:), intent(in) :: & system - character(len=*), intent(in) :: & - structure !< lattice structure + character(len=2), intent(in) :: & + lattice !< Bravais lattice (Pearson symbol) real(pReal), intent(in) :: & cOverA real(pReal), dimension(3,3,sum(active)) :: & @@ -1895,10 +1846,10 @@ function buildCoordinateSystem(active,potential,system,structure,cOverA) f, & !< index of my family s !< index of my system in current family - if (structure == 'tI' .and. cOverA > 2.0_pReal) & - call IO_error(131,ext_msg='buildCoordinateSystem:'//trim(structure)) - if (structure == 'hP' .and. (cOverA < 1.0_pReal .or. cOverA > 2.0_pReal)) & - call IO_error(131,ext_msg='buildCoordinateSystem:'//trim(structure)) + if (lattice == 'tI' .and. cOverA > 2.0_pReal) & + call IO_error(131,ext_msg='buildCoordinateSystem:'//trim(lattice)) + if (lattice == 'hP' .and. (cOverA < 1.0_pReal .or. cOverA > 2.0_pReal)) & + call IO_error(131,ext_msg='buildCoordinateSystem:'//trim(lattice)) a = 0 activeFamilies: do f = 1,size(active,1) @@ -1906,7 +1857,7 @@ function buildCoordinateSystem(active,potential,system,structure,cOverA) a = a + 1 p = sum(potential(1:f-1))+s - select case(structure) + select case(lattice) case ('cF','cI','tI') direction = system(1:3,p) @@ -1921,7 +1872,7 @@ function buildCoordinateSystem(active,potential,system,structure,cOverA) system(8,p)/cOverA ] ! plane (hkil)->(h (h+2k)/sqrt(3) l/(p/a)) case default - call IO_error(137,ext_msg='buildCoordinateSystem: '//trim(structure)) + call IO_error(137,ext_msg='buildCoordinateSystem: '//trim(lattice)) end select @@ -1950,9 +1901,9 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc) Q, & !< Total rotation: Q = R*B S !< Eigendeformation tensor for phase transformation real(pReal), intent(in) :: & - cOverA, & !< c/a for target hex structure - a_bcc, & !< lattice parameter a for target bcc structure - a_fcc !< lattice parameter a for parent fcc structure + cOverA, & !< c/a for target hex lattice + a_bcc, & !< lattice parameter a for bcc target lattice + a_fcc !< lattice parameter a for fcc parent lattice type(rotation) :: & R, & !< Pitsch rotation @@ -2126,9 +2077,10 @@ end function getlabels function lattice_equivalent_nu(C,assumption) result(nu) real(pReal), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation) - character(len=*), intent(in) :: assumption !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress) - real(pReal) :: K, mu, nu + character(len=5), intent(in) :: assumption !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress) + real(pReal) :: nu + real(pReal) :: K, mu logical :: error real(pReal), dimension(6,6) :: S @@ -2158,8 +2110,8 @@ end function lattice_equivalent_nu function lattice_equivalent_mu(C,assumption) result(mu) real(pReal), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation) - character(len=*), intent(in) :: assumption !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress) - real(pReal) :: mu + character(len=5), intent(in) :: assumption !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress) + real(pReal) :: mu logical :: error real(pReal), dimension(6,6) :: S @@ -2203,10 +2155,10 @@ subroutine selfTest do i = 1, 10 call random_number(C) - C_cF = applyLatticeSymmetryC66(C,'cI') - C_cI = applyLatticeSymmetryC66(C,'cF') - C_hP = applyLatticeSymmetryC66(C,'hP') - C_tI = applyLatticeSymmetryC66(C,'tI') + C_cF = lattice_symmetrize_C66(C,'cI') + C_cI = lattice_symmetrize_C66(C,'cF') + C_hP = lattice_symmetrize_C66(C,'hP') + C_tI = lattice_symmetrize_C66(C,'tI') if (any(dNeq(C_cI,transpose(C_cF)))) error stop 'SymmetryC66/cI-cF' if (any(dNeq(C_cF,transpose(C_cI)))) error stop 'SymmetryC66/cF-cI' @@ -2226,10 +2178,10 @@ subroutine selfTest if (any(dNeq(C(4,4),[C_tI(4,4),C_tI(5,5)]))) error stop 'SymmetryC_44-55/tI' call random_number(T) - T_cF = lattice_applyLatticeSymmetry33(T,'cI') - T_cI = lattice_applyLatticeSymmetry33(T,'cF') - T_hP = lattice_applyLatticeSymmetry33(T,'hP') - T_tI = lattice_applyLatticeSymmetry33(T,'tI') + T_cF = lattice_symmetrize_33(T,'cI') + T_cI = lattice_symmetrize_33(T,'cF') + T_hP = lattice_symmetrize_33(T,'hP') + T_tI = lattice_symmetrize_33(T,'tI') if (any(dNeq0(T_cF) .and. math_I3<1.0_pReal)) error stop 'Symmetry33/c' if (any(dNeq0(T_hP) .and. math_I3<1.0_pReal)) error stop 'Symmetry33/hP' @@ -2244,7 +2196,7 @@ subroutine selfTest call random_number(C) C(1,1) = C(1,1) + C(1,2) + 0.1_pReal C(4,4) = 0.5_pReal * (C(1,1) - C(1,2)) - C = applyLatticeSymmetryC66(C,'cI') + C = lattice_symmetrize_C66(C,'cI') if(dNeq(C(4,4),lattice_equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/voigt' if(dNeq(C(4,4),lattice_equivalent_mu(C,'reuss'),1.0e-12_pReal)) error stop 'equivalent_mu/reuss' diff --git a/src/math.f90 b/src/math.f90 index 191791ee2..36fe91bac 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -1068,6 +1068,7 @@ integer pure function math_factorial(n) integer, intent(in) :: n + math_factorial = product(math_range(n)) end function math_factorial @@ -1081,6 +1082,7 @@ integer pure function math_binomial(n,k) integer, intent(in) :: n, k integer :: i, k_, n_ + k_ = min(k,n-k) n_ = n math_binomial = merge(1,0,k_>-1) ! handling special cases k < 0 or k > n @@ -1095,15 +1097,13 @@ end function math_binomial !-------------------------------------------------------------------------------------------------- !> @brief multinomial coefficient !-------------------------------------------------------------------------------------------------- -integer pure function math_multinomial(alpha) +integer pure function math_multinomial(k) - integer, intent(in), dimension(:) :: alpha + integer, intent(in), dimension(:) :: k integer :: i - math_multinomial = 1 - do i = 1, size(alpha) - math_multinomial = math_multinomial*math_binomial(sum(alpha(1:i)),alpha(i)) - enddo + + math_multinomial = product([(math_binomial(sum(k(1:i)),k(i)),i=1,size(k))]) end function math_multinomial @@ -1116,6 +1116,7 @@ real(pReal) pure function math_volTetrahedron(v1,v2,v3,v4) real(pReal), dimension (3), intent(in) :: v1,v2,v3,v4 real(pReal), dimension (3,3) :: m + m(1:3,1) = v1-v2 m(1:3,2) = v1-v3 m(1:3,3) = v1-v4 @@ -1289,6 +1290,9 @@ subroutine selfTest if(math_binomial(49,6) /= 13983816) & error stop 'math_binomial' + if(math_multinomial([1,2,3,4]) /= 12600) & + error stop 'math_multinomial' + ijk = cshift([1,2,3],int(r*1.0e2_pReal)) if(dNeq(math_LeviCivita(ijk(1),ijk(2),ijk(3)),+1.0_pReal)) & error stop 'math_LeviCivita(even)' diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 6765d3d0d..5bb6d5f44 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -27,17 +27,17 @@ module FEM_utilities logical, public :: cutBack = .false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill integer, public, parameter :: maxFields = 6 integer, public :: nActiveFields = 0 - + !-------------------------------------------------------------------------------------------------- ! grid related information information real(pReal), public :: wgt !< weighting factor 1/Nelems - + !-------------------------------------------------------------------------------------------------- ! field labels information character(len=*), parameter, public :: & FIELD_MECH_label = 'mechanical' - + enum, bind(c); enumerator :: & FIELD_UNDEFINED_ID, & FIELD_MECH_ID @@ -48,7 +48,7 @@ module FEM_utilities COMPONENT_MECH_Y_ID, & COMPONENT_MECH_Z_ID end enum - + !-------------------------------------------------------------------------------------------------- ! variables controlling debugging logical :: & @@ -57,23 +57,23 @@ module FEM_utilities !-------------------------------------------------------------------------------------------------- ! derived types type, public :: tSolutionState !< return type of solution from FEM solver variants - logical :: converged = .true. - logical :: stagConverged = .true. + logical :: converged = .true. + logical :: stagConverged = .true. integer :: iterationsNeeded = 0 end type tSolutionState - + type, public :: tComponentBC integer(kind(COMPONENT_UNDEFINED_ID)) :: ID real(pReal), allocatable, dimension(:) :: Value - logical, allocatable, dimension(:) :: Mask + logical, allocatable, dimension(:) :: Mask end type tComponentBC - + type, public :: tFieldBC integer(kind(FIELD_UNDEFINED_ID)) :: ID integer :: nComponents = 0 type(tComponentBC), allocatable :: componentBC(:) end type tFieldBC - + type, public :: tLoadCase real(pReal) :: time = 0.0_pReal !< length of increment integer :: incs = 0, & !< number of increments @@ -83,7 +83,7 @@ module FEM_utilities integer, allocatable, dimension(:) :: faceID type(tFieldBC), allocatable, dimension(:) :: fieldBC end type tLoadCase - + public :: & FEM_utilities_init, & utilities_constitutiveResponse, & @@ -94,14 +94,14 @@ module FEM_utilities COMPONENT_MECH_Y_ID, & COMPONENT_MECH_Z_ID -contains +contains !ToDo: use functions in variable call !-------------------------------------------------------------------------------------------------- !> @brief allocates all neccessary fields, sets debug flags !-------------------------------------------------------------------------------------------------- subroutine FEM_utilities_init - + character(len=pStringLen) :: petsc_optionsOrder class(tNode), pointer :: & num_mesh, & @@ -113,7 +113,7 @@ subroutine FEM_utilities_init PetscErrorCode :: ierr print'(/,a)', ' <<<+- FEM_utilities init -+>>>' - + num_mesh => config_numerics%get('mesh',defaultVal=emptyDict) structOrder = num_mesh%get_asInt('structOrder', defaultVal = 2) @@ -141,7 +141,7 @@ subroutine FEM_utilities_init write(petsc_optionsOrder,'(a,i0)') '-mechFE_petscspace_degree ', structOrder call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsOrder),ierr) CHKERRQ(ierr) - + wgt = 1.0/real(mesh_maxNips*mesh_NcpElemsGlobal,pReal) @@ -152,20 +152,21 @@ end subroutine FEM_utilities_init !> @brief calculates constitutive response !-------------------------------------------------------------------------------------------------- subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData) - + real(pReal), intent(in) :: timeinc !< loading time logical, intent(in) :: forwardData !< age results - + real(pReal),intent(out), dimension(3,3) :: P_av !< average PK stress - + PetscErrorCode :: ierr print'(/,a)', ' ... evaluating constitutive response ......................................' - call materialpoint_stressAndItsTangent(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) ! calculate P field + call homogenization_mechanical_response(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) ! calculate P field + if (.not. terminallyIll) & + call homogenization_mechanical_response2(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) + cutBack = .false. - cutBack = .false. ! reset cutBack status - P_av = sum(homogenization_P,dim=3) * wgt call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr) @@ -198,12 +199,12 @@ subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCVa do dof = offset+comp+1, offset+numDof, numComp localArray(dof) = localArray(dof) + BCValue + BCDotValue*timeinc enddo - enddo + enddo call VecRestoreArrayF90(localVec,localArray,ierr); CHKERRQ(ierr) call VecAssemblyBegin(localVec, ierr); CHKERRQ(ierr) call VecAssemblyEnd (localVec, ierr); CHKERRQ(ierr) if (nBcPoints > 0) call ISRestoreIndicesF90(bcPointsIS,bcPoints,ierr) - + end subroutine utilities_projectBCValues end module FEM_utilities diff --git a/src/phase.f90 b/src/phase.f90 index e396e5661..d90f41ccc 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -108,9 +108,13 @@ module phase logical, intent(in) :: includeL end subroutine mechanical_restore + module subroutine damage_restore(ce) + integer, intent(in) :: ce + end subroutine damage_restore - module function phase_mechanical_dPdF(dt,co,ce) result(dPdF) - real(pReal), intent(in) :: dt + + module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF) + real(pReal), intent(in) :: Delta_t integer, intent(in) :: & co, & !< counter in constituent loop ce @@ -164,8 +168,8 @@ module phase real(pReal) :: dot_T end function thermal_dot_T - module function damage_phi(ph,me) result(phi) - integer, intent(in) :: ph,me + module function damage_phi(ph,en) result(phi) + integer, intent(in) :: ph,en real(pReal) :: phi end function damage_phi @@ -209,29 +213,27 @@ module phase ! == cleaned:end =================================================================================== - module function thermal_stress(Delta_t,ph,en) result(converged_) + module function phase_thermal_constitutive(Delta_t,ph,en) result(converged_) real(pReal), intent(in) :: Delta_t integer, intent(in) :: ph, en logical :: converged_ - end function thermal_stress + end function phase_thermal_constitutive - module function integrateDamageState(dt,co,ce) result(broken) - real(pReal), intent(in) :: dt - integer, intent(in) :: & - ce, & - co - logical :: broken - end function integrateDamageState - - module function crystallite_stress(dt,co,ip,el) result(converged_) - real(pReal), intent(in) :: dt + module function phase_damage_constitutive(Delta_t,co,ip,el) result(converged_) + real(pReal), intent(in) :: Delta_t integer, intent(in) :: co, ip, el logical :: converged_ - end function crystallite_stress + end function phase_damage_constitutive - !ToDo: Try to merge the all stiffness functions + module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged_) + real(pReal), intent(in) :: Delta_t + integer, intent(in) :: co, ip, el + logical :: converged_ + end function phase_mechanical_constitutive + + !ToDo: Merge all the stiffness functions module function phase_homogenizedC(ph,en) result(C) integer, intent(in) :: ph, en real(pReal), dimension(6,6) :: C @@ -271,8 +273,8 @@ module phase end subroutine plastic_dependentState - module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) - integer, intent(in) :: ph, me + module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,en) + integer, intent(in) :: ph, en real(pReal), intent(in), dimension(3,3) :: & S real(pReal), intent(out), dimension(3,3) :: & @@ -315,14 +317,14 @@ module phase plastic_nonlocal_updateCompatibility, & converged, & crystallite_init, & - crystallite_stress, & - thermal_stress, & + phase_mechanical_constitutive, & + phase_thermal_constitutive, & + phase_damage_constitutive, & phase_mechanical_dPdF, & crystallite_orientations, & crystallite_push33ToRef, & phase_restartWrite, & phase_restartRead, & - integrateDamageState, & phase_thermal_setField, & phase_set_phi, & phase_P, & @@ -435,17 +437,9 @@ subroutine phase_restore(ce,includeL) logical, intent(in) :: includeL integer, intent(in) :: ce - integer :: & - co - - - do co = 1,homogenization_Nconstituents(material_homogenizationID(ce)) - if (damageState(material_phaseID(co,ce))%sizeState > 0) & - damageState(material_phaseID(co,ce))%state( :,material_phaseEntry(co,ce)) = & - damageState(material_phaseID(co,ce))%state0(:,material_phaseEntry(co,ce)) - enddo call mechanical_restore(ce,includeL) + call damage_restore(ce) end subroutine phase_restore @@ -545,10 +539,6 @@ subroutine crystallite_init() phases => config_material%get('phase') - do ph = 1, phases%length - if (damageState(ph)%sizeState > 0) allocate(damageState(ph)%subState0,source=damageState(ph)%state0) ! ToDo: hack - enddo - print'(a42,1x,i10)', ' # of elements: ', eMax print'(a42,1x,i10)', ' # of integration points/element: ', iMax print'(a42,1x,i10)', 'max # of constituents/integration point: ', cMax diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index d8ac5046a..b69969898 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -5,7 +5,7 @@ submodule(phase) damage type :: tDamageParameters real(pReal) :: mu = 0.0_pReal !< viscosity - real(pReal), dimension(3,3) :: K = 0.0_pReal !< conductivity/diffusivity + real(pReal), dimension(3,3) :: D = 0.0_pReal !< conductivity/diffusivity end type tDamageParameters enum, bind(c); enumerator :: & @@ -18,7 +18,7 @@ submodule(phase) damage type :: tDataContainer - real(pReal), dimension(:), allocatable :: phi, d_phi_d_dot_phi + real(pReal), dimension(:), allocatable :: phi end type tDataContainer integer(kind(DAMAGE_UNDEFINED_ID)), dimension(:), allocatable :: & @@ -39,8 +39,8 @@ submodule(phase) damage end function isobrittle_init - module subroutine isobrittle_deltaState(C, Fe, ph, me) - integer, intent(in) :: ph,me + module subroutine isobrittle_deltaState(C, Fe, ph, en) + integer, intent(in) :: ph,en real(pReal), intent(in), dimension(3,3) :: & Fe real(pReal), intent(in), dimension(6,6) :: & @@ -48,8 +48,8 @@ submodule(phase) damage end subroutine isobrittle_deltaState - module subroutine anisobrittle_dotState(S, ph, me) - integer, intent(in) :: ph,me + module subroutine anisobrittle_dotState(S, ph, en) + integer, intent(in) :: ph,en real(pReal), intent(in), dimension(3,3) :: & S end subroutine anisobrittle_dotState @@ -97,7 +97,6 @@ module subroutine damage_init Nmembers = count(material_phaseID == ph) allocate(current(ph)%phi(Nmembers),source=1.0_pReal) - allocate(current(ph)%d_phi_d_dot_phi(Nmembers),source=0.0_pReal) phase => phases%get(ph) sources => phase%get('damage',defaultVal=emptyList) @@ -105,10 +104,10 @@ module subroutine damage_init if (sources%length == 1) then damage_active = .true. source => sources%get(1) - param(ph)%mu = source%get_asFloat('mu',defaultVal=0.0_pReal) ! ToDo: make mandatory? - param(ph)%K(1,1) = source%get_asFloat('K_11',defaultVal=0.0_pReal) ! ToDo: make mandatory? - param(ph)%K(3,3) = source%get_asFloat('K_33',defaultVal=0.0_pReal) ! ToDo: depends on symmetry - param(ph)%K = lattice_applyLatticeSymmetry33(param(ph)%K,phase_lattice(ph)) + param(ph)%mu = source%get_asFloat('mu') + param(ph)%D(1,1) = source%get_asFloat('D_11') + if (any(phase_lattice(ph) == ['hP','tI'])) param(ph)%D(3,3) = source%get_asFloat('D_33') + param(ph)%D = lattice_symmetrize_33(param(ph)%D,phase_lattice(ph)) endif enddo @@ -125,6 +124,29 @@ module subroutine damage_init end subroutine damage_init +!-------------------------------------------------------------------------------------------------- +!> @brief calculate stress (P) +!-------------------------------------------------------------------------------------------------- +module function phase_damage_constitutive(Delta_t,co,ip,el) result(converged_) + + real(pReal), intent(in) :: Delta_t + integer, intent(in) :: & + co, & + ip, & + el + logical :: converged_ + + integer :: & + ph, en + + ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) + en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) + + converged_ = .not. integrateDamageState(Delta_t,ph,en) + +end function phase_damage_constitutive + + !-------------------------------------------------------------------------------------------------- !> @brief returns the degraded/modified elasticity matrix !-------------------------------------------------------------------------------------------------- @@ -144,6 +166,26 @@ module function phase_damage_C(C_homogenized,ph,en) result(C) end function phase_damage_C +!-------------------------------------------------------------------------------------------------- +!> @brief Restore data after homog cutback. +!-------------------------------------------------------------------------------------------------- +module subroutine damage_restore(ce) + + integer, intent(in) :: ce + + integer :: & + co + + + do co = 1,homogenization_Nconstituents(material_homogenizationID(ce)) + if (damageState(material_phaseID(co,ce))%sizeState > 0) & + damageState(material_phaseID(co,ce))%state( :,material_phaseEntry(co,ce)) = & + damageState(material_phaseID(co,ce))%state0(:,material_phaseEntry(co,ce)) + enddo + +end subroutine damage_restore + + !---------------------------------------------------------------------------------------------- !< @brief returns local part of nonlocal damage driving force !---------------------------------------------------------------------------------------------- @@ -173,23 +215,20 @@ module function phase_f_phi(phi,co,ce) result(f) end function phase_f_phi - !-------------------------------------------------------------------------------------------------- !> @brief integrate stress, state with adaptive 1st order explicit Euler method !> using Fixed Point Iteration to adapt the stepsize !-------------------------------------------------------------------------------------------------- -module function integrateDamageState(dt,co,ce) result(broken) +function integrateDamageState(Delta_t,ph,en) result(broken) - real(pReal), intent(in) :: dt + real(pReal), intent(in) :: Delta_t integer, intent(in) :: & - ce, & - co + ph, & + en logical :: broken integer :: & NiterationState, & !< number of iterations in state loop - ph, & - me, & size_so real(pReal) :: & zeta @@ -199,8 +238,6 @@ module function integrateDamageState(dt,co,ce) result(broken) logical :: & converged_ - ph = material_phaseID(co,ce) - me = material_phaseEntry(co,ce) if (damageState(ph)%sizeState == 0) then broken = .false. @@ -208,37 +245,37 @@ module function integrateDamageState(dt,co,ce) result(broken) endif converged_ = .true. - broken = phase_damage_collectDotState(ph,me) + broken = phase_damage_collectDotState(ph,en) if(broken) return - size_so = damageState(ph)%sizeDotState - damageState(ph)%state(1:size_so,me) = damageState(ph)%subState0(1:size_so,me) & - + damageState(ph)%dotState (1:size_so,me) * dt - source_dotState(1:size_so,2) = 0.0_pReal + size_so = damageState(ph)%sizeDotState + damageState(ph)%state(1:size_so,en) = damageState(ph)%state0 (1:size_so,en) & + + damageState(ph)%dotState(1:size_so,en) * Delta_t + source_dotState(1:size_so,2) = 0.0_pReal iteration: do NiterationState = 1, num%nState if(nIterationState > 1) source_dotState(1:size_so,2) = source_dotState(1:size_so,1) - source_dotState(1:size_so,1) = damageState(ph)%dotState(:,me) + source_dotState(1:size_so,1) = damageState(ph)%dotState(:,en) - broken = phase_damage_collectDotState(ph,me) + broken = phase_damage_collectDotState(ph,en) if(broken) exit iteration - zeta = damper(damageState(ph)%dotState(:,me),source_dotState(1:size_so,1),source_dotState(1:size_so,2)) - damageState(ph)%dotState(:,me) = damageState(ph)%dotState(:,me) * zeta & + zeta = damper(damageState(ph)%dotState(:,en),source_dotState(1:size_so,1),source_dotState(1:size_so,2)) + damageState(ph)%dotState(:,en) = damageState(ph)%dotState(:,en) * zeta & + source_dotState(1:size_so,1)* (1.0_pReal - zeta) - r(1:size_so) = damageState(ph)%state (1:size_so,me) & - - damageState(ph)%subState0(1:size_so,me) & - - damageState(ph)%dotState (1:size_so,me) * dt - damageState(ph)%state(1:size_so,me) = damageState(ph)%state(1:size_so,me) - r(1:size_so) + r(1:size_so) = damageState(ph)%state (1:size_so,en) & + - damageState(ph)%State0 (1:size_so,en) & + - damageState(ph)%dotState(1:size_so,en) * Delta_t + damageState(ph)%state(1:size_so,en) = damageState(ph)%state(1:size_so,en) - r(1:size_so) converged_ = converged_ .and. converged(r(1:size_so), & - damageState(ph)%state(1:size_so,me), & + damageState(ph)%state(1:size_so,en), & damageState(ph)%atol(1:size_so)) if(converged_) then - broken = phase_damage_deltaState(mechanical_F_e(ph,me),ph,me) + broken = phase_damage_deltaState(mechanical_F_e(ph,en),ph,en) exit iteration endif @@ -300,11 +337,11 @@ end subroutine damage_results !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -function phase_damage_collectDotState(ph,me) result(broken) +function phase_damage_collectDotState(ph,en) result(broken) integer, intent(in) :: & ph, & - me !< counter in source loop + en !< counter in source loop logical :: broken @@ -315,11 +352,11 @@ function phase_damage_collectDotState(ph,me) result(broken) sourceType: select case (phase_damage(ph)) case (DAMAGE_ANISOBRITTLE_ID) sourceType - call anisobrittle_dotState(mechanical_S(ph,me), ph,me) ! correct stress? + call anisobrittle_dotState(mechanical_S(ph,en), ph,en) ! correct stress? end select sourceType - broken = broken .or. any(IEEE_is_NaN(damageState(ph)%dotState(:,me))) + broken = broken .or. any(IEEE_is_NaN(damageState(ph)%dotState(:,en))) endif @@ -347,9 +384,10 @@ module function phase_K_phi(co,ce) result(K) integer, intent(in) :: co, ce real(pReal), dimension(3,3) :: K - - - K = crystallite_push33ToRef(co,ce,param(material_phaseID(co,ce))%K) + real(pReal), parameter :: l = 1.0_pReal + + K = crystallite_push33ToRef(co,ce,param(material_phaseID(co,ce))%D) \ + * l**2.0_pReal end function phase_K_phi @@ -358,11 +396,11 @@ end function phase_K_phi !> @brief for constitutive models having an instantaneous change of state !> will return false if delta state is not needed/supported by the constitutive model !-------------------------------------------------------------------------------------------------- -function phase_damage_deltaState(Fe, ph, me) result(broken) +function phase_damage_deltaState(Fe, ph, en) result(broken) integer, intent(in) :: & ph, & - me + en real(pReal), intent(in), dimension(3,3) :: & Fe !< elastic deformation gradient integer :: & @@ -379,13 +417,13 @@ function phase_damage_deltaState(Fe, ph, me) result(broken) sourceType: select case (phase_damage(ph)) case (DAMAGE_ISOBRITTLE_ID) sourceType - call isobrittle_deltaState(phase_homogenizedC(ph,me), Fe, ph,me) - broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,me))) + call isobrittle_deltaState(phase_homogenizedC(ph,en), Fe, ph,en) + broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,en))) if(.not. broken) then myOffset = damageState(ph)%offsetDeltaState mySize = damageState(ph)%sizeDeltaState - damageState(ph)%state(myOffset + 1: myOffset + mySize,me) = & - damageState(ph)%state(myOffset + 1: myOffset + mySize,me) + damageState(ph)%deltaState(1:mySize,me) + damageState(ph)%state(myOffset + 1: myOffset + mySize,en) = & + damageState(ph)%state(myOffset + 1: myOffset + mySize,en) + damageState(ph)%deltaState(1:mySize,en) endif end select sourceType @@ -436,13 +474,13 @@ module subroutine phase_set_phi(phi,co,ce) end subroutine phase_set_phi -module function damage_phi(ph,me) result(phi) +module function damage_phi(ph,en) result(phi) - integer, intent(in) :: ph, me + integer, intent(in) :: ph, en real(pReal) :: phi - phi = current(ph)%phi(me) + phi = current(ph)%phi(en) end function damage_phi diff --git a/src/phase_damage_anisobrittle.f90 b/src/phase_damage_anisobrittle.f90 index 9984e5514..e8f4e7f98 100644 --- a/src/phase_damage_anisobrittle.f90 +++ b/src/phase_damage_anisobrittle.f90 @@ -40,7 +40,7 @@ module function anisobrittle_init() result(mySources) phase, & sources, & src - integer :: Nmembers,p + integer :: Nmembers,ph integer, dimension(:), allocatable :: N_cl character(len=pStringLen) :: extmsg = '' @@ -56,12 +56,12 @@ module function anisobrittle_init() result(mySources) allocate(param(phases%length)) - do p = 1, phases%length - if(mySources(p)) then - phase => phases%get(p) - sources => phase%get('damage') + do ph = 1, phases%length + if(mySources(ph)) then + phase => phases%get(ph) + sources => phase%get('damage') - associate(prm => param(p)) + associate(prm => param(ph)) src => sources%get(1) N_cl = src%get_as1dInt('N_cl',defaultVal=emptyIntArray) @@ -73,8 +73,7 @@ module function anisobrittle_init() result(mySources) prm%s_crit = src%get_as1dFloat('s_crit', requiredSize=size(N_cl)) prm%g_crit = src%get_as1dFloat('g_crit', requiredSize=size(N_cl)) - prm%cleavage_systems = lattice_SchmidMatrix_cleavage(N_cl,phase%get_asString('lattice'),& - phase%get_asFloat('c/a',defaultVal=0.0_pReal)) + prm%cleavage_systems = lattice_SchmidMatrix_cleavage(N_cl,phase_lattice(ph),phase_cOverA(ph)) ! expand: family => system prm%s_crit = math_expand(prm%s_crit,N_cl) @@ -92,17 +91,15 @@ module function anisobrittle_init() result(mySources) if (any(prm%g_crit < 0.0_pReal)) extmsg = trim(extmsg)//' g_crit' if (any(prm%s_crit < 0.0_pReal)) extmsg = trim(extmsg)//' s_crit' - Nmembers = count(material_phaseID==p) - call phase_allocateState(damageState(p),Nmembers,1,1,0) - damageState(p)%atol = src%get_asFloat('atol_phi',defaultVal=1.0e-9_pReal) - if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' atol_phi' + Nmembers = count(material_phaseID==ph) + call phase_allocateState(damageState(ph),Nmembers,1,1,0) + damageState(ph)%atol = src%get_asFloat('atol_phi',defaultVal=1.0e-9_pReal) + if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' atol_phi' - end associate + end associate -!-------------------------------------------------------------------------------------------------- -! exit if any parameter is out of range - if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_anisoBrittle)') - endif + if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_anisoBrittle)') + endif enddo @@ -112,10 +109,10 @@ end function anisobrittle_init !-------------------------------------------------------------------------------------------------- !> @brief !-------------------------------------------------------------------------------------------------- -module subroutine anisobrittle_dotState(S, ph,me) +module subroutine anisobrittle_dotState(S, ph,en) integer, intent(in) :: & - ph,me + ph,en real(pReal), intent(in), dimension(3,3) :: & S @@ -126,15 +123,15 @@ module subroutine anisobrittle_dotState(S, ph,me) associate(prm => param(ph)) - damageState(ph)%dotState(1,me) = 0.0_pReal + damageState(ph)%dotState(1,en) = 0.0_pReal do i = 1, prm%sum_N_cl traction_d = math_tensordot(S,prm%cleavage_systems(1:3,1:3,1,i)) traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i)) traction_n = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i)) - traction_crit = prm%g_crit(i)*damage_phi(ph,me)**2.0_pReal + traction_crit = prm%g_crit(i)*damage_phi(ph,en)**2.0_pReal - damageState(ph)%dotState(1,me) = damageState(ph)%dotState(1,me) & + damageState(ph)%dotState(1,en) = damageState(ph)%dotState(1,en) & + prm%dot_o / prm%s_crit(i) & * ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**prm%q + & (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**prm%q + & @@ -171,10 +168,10 @@ end subroutine anisobrittle_results !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) +module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,en) integer, intent(in) :: & - ph,me + ph,en real(pReal), intent(in), dimension(3,3) :: & S real(pReal), intent(out), dimension(3,3) :: & @@ -193,7 +190,7 @@ module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) dLd_dTstar = 0.0_pReal associate(prm => param(ph)) do i = 1,prm%sum_N_cl - traction_crit = prm%g_crit(i)*damage_phi(ph,me)**2.0_pReal + traction_crit = prm%g_crit(i)*damage_phi(ph,en)**2.0_pReal traction_d = math_tensordot(S,prm%cleavage_systems(1:3,1:3,1,i)) if (abs(traction_d) > traction_crit + tol_math_check) then diff --git a/src/phase_damage_isobrittle.f90 b/src/phase_damage_isobrittle.f90 index aed6a91e9..066cab47e 100644 --- a/src/phase_damage_isobrittle.f90 +++ b/src/phase_damage_isobrittle.f90 @@ -13,7 +13,15 @@ submodule(phase:damage) isobrittle output end type tParameters - type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances) + type :: tIsobrittleState + real(pReal), pointer, dimension(:) :: & !< vectors along Nmembers + r_W !< ratio between actual and critical strain energy density + end type tIsobrittleState + + type(tParameters), allocatable, dimension(:) :: param !< containers of constitutive parameters (len Ninstances) + type(tIsobrittleState), allocatable, dimension(:) :: & + deltaState, & + state contains @@ -44,13 +52,15 @@ module function isobrittle_init() result(mySources) phases => config_material%get('phase') allocate(param(phases%length)) + allocate(state(phases%length)) + allocate(deltaState(phases%length)) do ph = 1, phases%length if(mySources(ph)) then - phase => phases%get(ph) - sources => phase%get('damage') + phase => phases%get(ph) + sources => phase%get('damage') - associate(prm => param(ph)) + associate(prm => param(ph), dlt => deltaState(ph), stt => state(ph)) src => sources%get(1) prm%W_crit = src%get_asFloat('W_crit') @@ -69,12 +79,14 @@ module function isobrittle_init() result(mySources) damageState(ph)%atol = src%get_asFloat('atol_phi',defaultVal=1.0e-9_pReal) if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' atol_phi' - end associate + stt%r_W => damageState(ph)%state(1,:) + dlt%r_W => damageState(ph)%deltaState(1,:) -!-------------------------------------------------------------------------------------------------- -! exit if any parameter is out of range - if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_isobrittle)') - endif + end associate + + + if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_isobrittle)') + endif enddo @@ -85,29 +97,27 @@ end function isobrittle_init !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -module subroutine isobrittle_deltaState(C, Fe, ph,me) +module subroutine isobrittle_deltaState(C, Fe, ph,en) - integer, intent(in) :: ph,me + integer, intent(in) :: ph,en real(pReal), intent(in), dimension(3,3) :: & Fe real(pReal), intent(in), dimension(6,6) :: & C real(pReal), dimension(6) :: & - strain + epsilon real(pReal) :: & - strainenergy + r_W - strain = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3) + epsilon = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3) - associate(prm => param(ph)) - strainenergy = 2.0_pReal*sum(strain*matmul(C,strain))/prm%W_crit - ! ToDo: check strainenergy = 2.0_pReal*dot_product(strain,matmul(C,strain))/prm%W_crit + associate(prm => param(ph), stt => state(ph), dlt => deltaState(ph)) + + r_W = 2.0_pReal*dot_product(epsilon,matmul(C,epsilon))/prm%W_crit + dlt%r_W(en) = merge(r_W - stt%r_W(en), 0.0_pReal, r_W > stt%r_W(en)) - damageState(ph)%deltaState(1,me) = merge(strainenergy - damageState(ph)%state(1,me), & - damageState(ph)%subState0(1,me) - damageState(ph)%state(1,me), & - strainenergy > damageState(ph)%subState0(1,me)) end associate end subroutine isobrittle_deltaState @@ -124,14 +134,15 @@ module subroutine isobrittle_results(phase,group) integer :: o - associate(prm => param(phase), & - stt => damageState(phase)%state) + associate(prm => param(phase), stt => damageState(phase)%state) ! point to state and output r_W (is scalar, not 1D vector) + outputsLoop: do o = 1,size(prm%output) select case(trim(prm%output(o))) case ('f_phi') - call results_writeDataset(stt,group,trim(prm%output(o)),'driving force','J/m³') + call results_writeDataset(stt,group,trim(prm%output(o)),'driving force','J/m³') ! Wrong, this is dimensionless end select enddo outputsLoop + end associate end subroutine isobrittle_results diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 7562c055f..c972ed3cd 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -44,9 +44,9 @@ submodule(phase) mechanical interface - module subroutine eigendeformation_init(phases) + module subroutine eigen_init(phases) class(tNode), pointer :: phases - end subroutine eigendeformation_init + end subroutine eigen_init module subroutine elastic_init(phases) class(tNode), pointer :: phases @@ -168,6 +168,20 @@ submodule(phase) mechanical integer, intent(in) :: ph,en end function plastic_dislotwin_homogenizedC + module function elastic_C66(ph) result(C66) + real(pReal), dimension(6,6) :: C66 + integer, intent(in) :: ph + end function elastic_C66 + + module function elastic_mu(ph) result(mu) + real(pReal) :: mu + integer, intent(in) :: ph + end function elastic_mu + + module function elastic_nu(ph) result(nu) + real(pReal) :: nu + integer, intent(in) :: ph + end function elastic_nu end interface type :: tOutput !< new requested output (per phase) @@ -197,8 +211,7 @@ module subroutine mechanical_init(materials,phases) Nmembers class(tNode), pointer :: & num_crystallite, & - material, & - constituents, & + phase, & mech @@ -303,7 +316,7 @@ module subroutine mechanical_init(materials,phases) end select - call eigendeformation_init(phases) + call eigen_init(phases) end subroutine mechanical_init @@ -977,9 +990,9 @@ end subroutine mechanical_forward !-------------------------------------------------------------------------------------------------- !> @brief calculate stress (P) !-------------------------------------------------------------------------------------------------- -module function crystallite_stress(dt,co,ip,el) result(converged_) +module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged_) - real(pReal), intent(in) :: dt + real(pReal), intent(in) :: Delta_t integer, intent(in) :: & co, & ip, & @@ -1009,17 +1022,13 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) subLi0 = phase_mechanical_Li0(ph)%data(1:3,1:3,en) subLp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,en) subState0 = plasticState(ph)%State0(:,en) - - if (damageState(ph)%sizeState > 0) & - damageState(ph)%subState0(:,en) = damageState(ph)%state0(:,en) - subFp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,en) subFi0 = phase_mechanical_Fi0(ph)%data(1:3,1:3,en) subF0 = phase_mechanical_F0(ph)%data(1:3,1:3,en) subFrac = 0.0_pReal - subStep = 1.0_pReal/num%subStepSizeCryst todo = .true. - converged_ = .false. ! pretend failed step of 1/subStepSizeCryst + subStep = 1.0_pReal/num%subStepSizeCryst + converged_ = .false. ! pretend failed step of 1/subStepSizeCryst todo = .true. cutbackLooping: do while (todo) @@ -1038,9 +1047,6 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) subFp0 = phase_mechanical_Fp(ph)%data(1:3,1:3,en) subFi0 = phase_mechanical_Fi(ph)%data(1:3,1:3,en) subState0 = plasticState(ph)%state(:,en) - if (damageState(ph)%sizeState > 0) & - damageState(ph)%subState0(:,en) = damageState(ph)%state(:,en) - endif !-------------------------------------------------------------------------------------------------- ! cut back (reduced time and restore) @@ -1048,16 +1054,13 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) subStep = num%subStepSizeCryst * subStep phase_mechanical_Fp(ph)%data(1:3,1:3,en) = subFp0 phase_mechanical_Fi(ph)%data(1:3,1:3,en) = subFi0 - phase_mechanical_S(ph)%data(1:3,1:3,en) = phase_mechanical_S0(ph)%data(1:3,1:3,en) ! why no subS0 ? is S0 of any use? + phase_mechanical_S(ph)%data(1:3,1:3,en) = phase_mechanical_S0(ph)%data(1:3,1:3,en) ! why no subS0 ? is S0 of any use? if (subStep < 1.0_pReal) then ! actual (not initial) cutback phase_mechanical_Lp(ph)%data(1:3,1:3,en) = subLp0 phase_mechanical_Li(ph)%data(1:3,1:3,en) = subLi0 endif plasticState(ph)%state(:,en) = subState0 - if (damageState(ph)%sizeState > 0) & - damageState(ph)%state(:,en) = damageState(ph)%subState0(:,en) - - todo = subStep > num%subStepMinCryst ! still on track or already done (beyond repair) + todo = subStep > num%subStepMinCryst ! still on track or already done (beyond repair) endif !-------------------------------------------------------------------------------------------------- @@ -1065,13 +1068,12 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) if (todo) then subF = subF0 & + subStep * (phase_mechanical_F(ph)%data(1:3,1:3,en) - phase_mechanical_F0(ph)%data(1:3,1:3,en)) - converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * dt,co,ip,el) - converged_ = converged_ .and. .not. integrateDamageState(subStep * dt,co,(el-1)*discretization_nIPs + ip) + converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * Delta_t,co,ip,el) endif enddo cutbackLooping -end function crystallite_stress +end function phase_mechanical_constitutive !-------------------------------------------------------------------------------------------------- @@ -1104,12 +1106,13 @@ module subroutine mechanical_restore(ce,includeL) end subroutine mechanical_restore + !-------------------------------------------------------------------------------------------------- !> @brief Calculate tangent (dPdF). !-------------------------------------------------------------------------------------------------- -module function phase_mechanical_dPdF(dt,co,ce) result(dPdF) +module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF) - real(pReal), intent(in) :: dt + real(pReal), intent(in) :: Delta_t integer, intent(in) :: & co, & !< counter in constituent loop ce @@ -1159,11 +1162,11 @@ module function phase_mechanical_dPdF(dt,co,ce) result(dPdF) lhs_3333 = 0.0_pReal; rhs_3333 = 0.0_pReal do o=1,3; do p=1,3 lhs_3333(1:3,1:3,o,p) = lhs_3333(1:3,1:3,o,p) & - + matmul(invSubFi0,dLidFi(1:3,1:3,o,p)) * dt + + matmul(invSubFi0,dLidFi(1:3,1:3,o,p)) * Delta_t lhs_3333(1:3,o,1:3,p) = lhs_3333(1:3,o,1:3,p) & + invFi*invFi(p,o) rhs_3333(1:3,1:3,o,p) = rhs_3333(1:3,1:3,o,p) & - - matmul(invSubFi0,dLidS(1:3,1:3,o,p)) * dt + - matmul(invSubFi0,dLidS(1:3,1:3,o,p)) * Delta_t enddo; enddo call math_invert(temp_99,error,math_3333to99(lhs_3333)) if (error) then @@ -1192,7 +1195,7 @@ module function phase_mechanical_dPdF(dt,co,ce) result(dPdF) temp_3333(1:3,1:3,p,o) = matmul(matmul(temp_33_2,dLpdS(1:3,1:3,p,o)), invFi) & + matmul(temp_33_3,dLidS(1:3,1:3,p,o)) enddo; enddo - lhs_3333 = math_mul3333xx3333(dSdFe,temp_3333) * dt & + lhs_3333 = math_mul3333xx3333(dSdFe,temp_3333) * Delta_t & + math_mul3333xx3333(dSdFi,dFidS) call math_invert(temp_99,error,math_eye(9)+math_3333to99(lhs_3333)) @@ -1208,7 +1211,7 @@ module function phase_mechanical_dPdF(dt,co,ce) result(dPdF) ! calculate dFpinvdF temp_3333 = math_mul3333xx3333(dLpdS,dSdF) do o=1,3; do p=1,3 - dFpinvdF(1:3,1:3,p,o) = - matmul(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi)) * dt + dFpinvdF(1:3,1:3,p,o) = - matmul(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi)) * Delta_t enddo; enddo !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_mechanical_eigen.f90 b/src/phase_mechanical_eigen.f90 index 019838689..d33adc5d1 100644 --- a/src/phase_mechanical_eigen.f90 +++ b/src/phase_mechanical_eigen.f90 @@ -32,7 +32,7 @@ submodule(phase:mechanical) eigen contains -module subroutine eigendeformation_init(phases) +module subroutine eigen_init(phases) class(tNode), pointer :: & phases @@ -68,7 +68,7 @@ module subroutine eigendeformation_init(phases) where(damage_anisobrittle_init()) model_damage = KINEMATICS_cleavage_opening_ID -end subroutine eigendeformation_init +end subroutine eigen_init !-------------------------------------------------------------------------------------------------- @@ -86,17 +86,17 @@ function kinematics_active(kinematics_label,kinematics_length) result(active_ki kinematics, & kinematics_type, & mechanics - integer :: p,k + integer :: ph,k phases => config_material%get('phase') allocate(active_kinematics(kinematics_length,phases%length), source = .false. ) - do p = 1, phases%length - phase => phases%get(p) + do ph = 1, phases%length + phase => phases%get(ph) mechanics => phase%get('mechanical') kinematics => mechanics%get('eigen',defaultVal=emptyList) do k = 1, kinematics%length kinematics_type => kinematics%get(k) - active_kinematics(k,p) = kinematics_type%get_asString('type') == kinematics_label + active_kinematics(k,ph) = kinematics_type%get_asString('type') == kinematics_label enddo enddo @@ -118,17 +118,17 @@ function kinematics_active2(kinematics_label) result(active_kinematics) phase, & kinematics, & kinematics_type - integer :: p + integer :: ph phases => config_material%get('phase') - allocate(active_kinematics(phases%length), source = .false. ) - do p = 1, phases%length - phase => phases%get(p) + allocate(active_kinematics(phases%length), source = .false.) + do ph = 1, phases%length + phase => phases%get(ph) kinematics => phase%get('damage',defaultVal=emptyList) if(kinematics%length < 1) return kinematics_type => kinematics%get(1) if (.not. kinematics_type%contains('type')) continue - active_kinematics(p) = kinematics_type%get_asString('type',defaultVal='n/a') == kinematics_label + active_kinematics(ph) = kinematics_type%get_asString('type',defaultVal='n/a') == kinematics_label enddo diff --git a/src/phase_mechanical_eigen_thermalexpansion.f90 b/src/phase_mechanical_eigen_thermalexpansion.f90 index dba02d70e..0976c1dfa 100644 --- a/src/phase_mechanical_eigen_thermalexpansion.f90 +++ b/src/phase_mechanical_eigen_thermalexpansion.f90 @@ -69,7 +69,7 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics) prm%A(3,3,3) = kinematic_type%get_asFloat('A_33,T^2',defaultVal=0.0_pReal) endif do i=1, size(prm%A,3) - prm%A(1:3,1:3,i) = lattice_applyLatticeSymmetry33(prm%A(1:3,1:3,i),phase_lattice(p)) + prm%A(1:3,1:3,i) = lattice_symmetrize_33(prm%A(1:3,1:3,i),phase_lattice(p)) enddo end associate endif diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 index 3249b2921..326c23c7a 100644 --- a/src/phase_mechanical_elastic.f90 +++ b/src/phase_mechanical_elastic.f90 @@ -2,7 +2,10 @@ submodule(phase:mechanical) elastic type :: tParameters real(pReal), dimension(6,6) :: & - C66 !< Elastic constants in Voig notation + C66 = 0.0_pReal !< Elastic constants in Voigt notation + real(pReal) :: & + mu, & + nu end type tParameters type(tParameters), allocatable, dimension(:) :: param @@ -37,6 +40,7 @@ module subroutine elastic_init(phases) if (elastic%get_asString('type') /= 'Hooke') call IO_error(200,ext_msg=elastic%get_asString('type')) associate(prm => param(ph)) + prm%C66(1,1) = elastic%get_asFloat('C_11') prm%C66(1,2) = elastic%get_asFloat('C_12') prm%C66(4,4) = elastic%get_asFloat('C_44') @@ -46,6 +50,14 @@ module subroutine elastic_init(phases) prm%C66(3,3) = elastic%get_asFloat('C_33') endif if (phase_lattice(ph) == 'tI') prm%C66(6,6) = elastic%get_asFloat('C_66') + + prm%C66 = lattice_symmetrize_C66(prm%C66,phase_lattice(ph)) + + prm%nu = lattice_equivalent_nu(prm%C66,'voigt') + prm%mu = lattice_equivalent_mu(prm%C66,'voigt') + + prm%C66 = math_sym3333to66(math_Voigt66to3333(prm%C66)) ! Literature data is in Voigt notation + end associate enddo @@ -104,10 +116,38 @@ module function phase_homogenizedC(ph,en) result(C) case (PLASTICITY_DISLOTWIN_ID) plasticType C = plastic_dislotwin_homogenizedC(ph,en) case default plasticType - C = lattice_C66(1:6,1:6,ph) + C = param(ph)%C66 end select plasticType end function phase_homogenizedC +module function elastic_C66(ph) result(C66) + real(pReal), dimension(6,6) :: C66 + integer, intent(in) :: ph + + + C66 = param(ph)%C66 + +end function elastic_C66 + +module function elastic_mu(ph) result(mu) + + real(pReal) :: mu + integer, intent(in) :: ph + + + mu = param(ph)%mu + +end function elastic_mu + +module function elastic_nu(ph) result(nu) + + real(pReal) :: nu + integer, intent(in) :: ph + + + nu = param(ph)%mu + +end function elastic_nu end submodule elastic diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index 0fa742b01..6ca8e440e 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -36,8 +36,8 @@ submodule(phase:plastic) dislotungsten forestProjection real(pReal), allocatable, dimension(:,:,:) :: & P_sl, & - nonSchmid_pos, & - nonSchmid_neg + P_nS_pos, & + P_nS_neg integer :: & sum_N_sl !< total number of active slip system character(len=pStringLen), allocatable, dimension(:) :: & @@ -129,30 +129,28 @@ module function plastic_dislotungsten_init() result(myPlasticity) prm%output = pl%get_as1dString('output',defaultVal=emptyStringArray) #endif - ! This data is read in already in lattice - prm%mu = lattice_mu(ph) + prm%mu = elastic_mu(ph) !-------------------------------------------------------------------------------------------------- ! slip related parameters N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray) prm%sum_N_sl = sum(abs(N_sl)) slipActive: if (prm%sum_N_sl > 0) then - prm%P_sl = lattice_SchmidMatrix_slip(N_sl,phase%get_asString('lattice'),& - phase%get_asFloat('c/a',defaultVal=0.0_pReal)) + prm%P_sl = lattice_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph)) - if(trim(phase%get_asString('lattice')) == 'cI') then + if (phase_lattice(ph) == 'cI') then a = pl%get_as1dFloat('a_nonSchmid',defaultVal = emptyRealArray) - prm%nonSchmid_pos = lattice_nonSchmidMatrix(N_sl,a,+1) - prm%nonSchmid_neg = lattice_nonSchmidMatrix(N_sl,a,-1) + prm%P_nS_pos = lattice_nonSchmidMatrix(N_sl,a,+1) + prm%P_nS_neg = lattice_nonSchmidMatrix(N_sl,a,-1) else - prm%nonSchmid_pos = prm%P_sl - prm%nonSchmid_neg = prm%P_sl + prm%P_nS_pos = prm%P_sl + prm%P_nS_neg = prm%P_sl endif prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dFloat('h_sl-sl'), & - phase%get_asString('lattice')) - prm%forestProjection = lattice_forestProjection_edge(N_sl,phase%get_asString('lattice'),& - phase%get_asFloat('c/a',defaultVal=0.0_pReal)) + phase_lattice(ph)) + prm%forestProjection = lattice_forestProjection_edge(N_sl,phase_lattice(ph),& + phase_cOverA(ph)) prm%forestProjection = transpose(prm%forestProjection) rho_mob_0 = pl%get_as1dFloat('rho_mob_0', requiredSize=size(N_sl)) @@ -163,39 +161,37 @@ module function plastic_dislotungsten_init() result(myPlasticity) prm%i_sl = pl%get_as1dFloat('i_sl', requiredSize=size(N_sl)) prm%tau_Peierls = pl%get_as1dFloat('tau_Peierls', requiredSize=size(N_sl)) - prm%p = pl%get_as1dFloat('p_sl', requiredSize=size(N_sl), & - defaultVal=[(1.0_pReal,i=1,size(N_sl))]) - prm%q = pl%get_as1dFloat('q_sl', requiredSize=size(N_sl), & - defaultVal=[(1.0_pReal,i=1,size(N_sl))]) + prm%p = pl%get_as1dFloat('p_sl', requiredSize=size(N_sl)) + prm%q = pl%get_as1dFloat('q_sl', requiredSize=size(N_sl)) prm%h = pl%get_as1dFloat('h', requiredSize=size(N_sl)) prm%w = pl%get_as1dFloat('w', requiredSize=size(N_sl)) prm%omega = pl%get_as1dFloat('omega', requiredSize=size(N_sl)) prm%B = pl%get_as1dFloat('B', requiredSize=size(N_sl)) - prm%D = pl%get_asFloat('D') - prm%D_0 = pl%get_asFloat('D_0') - prm%Q_cl = pl%get_asFloat('Q_cl') - prm%f_at = pl%get_asFloat('f_at') * prm%b_sl**3.0_pReal - prm%D_a = pl%get_asFloat('D_a') * prm%b_sl + prm%D = pl%get_asFloat('D') + prm%D_0 = pl%get_asFloat('D_0') + prm%Q_cl = pl%get_asFloat('Q_cl') + prm%f_at = pl%get_asFloat('f_at') * prm%b_sl**3.0_pReal + prm%D_a = pl%get_asFloat('D_a') * prm%b_sl prm%dipoleformation = .not. pl%get_asBool('no_dipole_formation', defaultVal = .false.) ! expand: family => system - rho_mob_0 = math_expand(rho_mob_0, N_sl) - rho_dip_0 = math_expand(rho_dip_0, N_sl) - prm%q = math_expand(prm%q, N_sl) - prm%p = math_expand(prm%p, N_sl) - prm%Q_s = math_expand(prm%Q_s, N_sl) - prm%b_sl = math_expand(prm%b_sl, N_sl) - prm%h = math_expand(prm%h, N_sl) - prm%w = math_expand(prm%w, N_sl) - prm%omega = math_expand(prm%omega, N_sl) - prm%tau_Peierls = math_expand(prm%tau_Peierls, N_sl) - prm%v_0 = math_expand(prm%v_0, N_sl) - prm%B = math_expand(prm%B, N_sl) - prm%i_sl = math_expand(prm%i_sl, N_sl) - prm%f_at = math_expand(prm%f_at, N_sl) - prm%D_a = math_expand(prm%D_a, N_sl) + rho_mob_0 = math_expand(rho_mob_0, N_sl) + rho_dip_0 = math_expand(rho_dip_0, N_sl) + prm%q = math_expand(prm%q, N_sl) + prm%p = math_expand(prm%p, N_sl) + prm%Q_s = math_expand(prm%Q_s, N_sl) + prm%b_sl = math_expand(prm%b_sl, N_sl) + prm%h = math_expand(prm%h, N_sl) + prm%w = math_expand(prm%w, N_sl) + prm%omega = math_expand(prm%omega, N_sl) + prm%tau_Peierls = math_expand(prm%tau_Peierls, N_sl) + prm%v_0 = math_expand(prm%v_0, N_sl) + prm%B = math_expand(prm%B, N_sl) + prm%i_sl = math_expand(prm%i_sl, N_sl) + prm%f_at = math_expand(prm%f_at, N_sl) + prm%D_a = math_expand(prm%D_a, N_sl) ! sanity checks if ( prm%D_0 <= 0.0_pReal) extmsg = trim(extmsg)//' D_0' @@ -298,8 +294,8 @@ pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp, & Lp = Lp + (dot_gamma_pos(i)+dot_gamma_neg(i))*prm%P_sl(1:3,1:3,i) forall (k=1:3,l=1:3,m=1:3,n=1:3) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & - + ddot_gamma_dtau_pos(i) * prm%P_sl(k,l,i) * prm%nonSchmid_pos(m,n,i) & - + ddot_gamma_dtau_neg(i) * prm%P_sl(k,l,i) * prm%nonSchmid_neg(m,n,i) + + ddot_gamma_dtau_pos(i) * prm%P_sl(k,l,i) * prm%P_nS_pos(m,n,i) & + + ddot_gamma_dtau_neg(i) * prm%P_sl(k,l,i) * prm%P_nS_neg(m,n,i) enddo end associate @@ -323,7 +319,7 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en) real(pReal) :: & VacancyDiffusion real(pReal), dimension(param(ph)%sum_N_sl) :: & - gdot_pos, gdot_neg,& + dot_gamma_pos, dot_gamma_neg,& tau_pos,& tau_neg, & v_cl, & @@ -335,10 +331,10 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en) dot => dotState(ph), dst => dependentState(ph)) call kinetics(Mp,T,ph,en,& - gdot_pos,gdot_neg, & + dot_gamma_pos,dot_gamma_neg, & tau_pos_out = tau_pos,tau_neg_out = tau_neg) - dot%gamma_sl(:,en) = (gdot_pos+gdot_neg) ! ToDo: needs to be abs + dot%gamma_sl(:,en) = (dot_gamma_pos+dot_gamma_neg) ! ToDo: needs to be abs VacancyDiffusion = prm%D_0*exp(-prm%Q_cl/(kB*T)) where(dEq0(tau_pos)) ! ToDo: use avg of pos and neg @@ -380,13 +376,14 @@ module subroutine dislotungsten_dependentState(ph,en) real(pReal), dimension(param(ph)%sum_N_sl) :: & dislocationSpacing + associate(prm => param(ph), stt => state(ph),dst => dependentState(ph)) - dislocationSpacing = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en))) - dst%threshold_stress(:,en) = prm%mu*prm%b_sl & - * sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en))) + dislocationSpacing = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en))) + dst%threshold_stress(:,en) = prm%mu*prm%b_sl & + * sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en))) - dst%Lambda_sl(:,en) = prm%D/(1.0_pReal+prm%D*dislocationSpacing/prm%i_sl) + dst%Lambda_sl(:,en) = prm%D/(1.0_pReal+prm%D*dislocationSpacing/prm%i_sl) end associate @@ -466,8 +463,8 @@ pure subroutine kinetics(Mp,T,ph,en, & associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) do j = 1, prm%sum_N_sl - tau_pos(j) = math_tensordot(Mp,prm%nonSchmid_pos(1:3,1:3,j)) - tau_neg(j) = math_tensordot(Mp,prm%nonSchmid_neg(1:3,1:3,j)) + tau_pos(j) = math_tensordot(Mp,prm%P_nS_pos(1:3,1:3,j)) + tau_neg(j) = math_tensordot(Mp,prm%P_nS_neg(1:3,1:3,j)) enddo diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index de0d20b11..6ecf9ee7a 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -184,27 +184,23 @@ module function plastic_dislotwin_init() result(myPlasticity) #endif ! This data is read in already in lattice - prm%mu = lattice_mu(ph) - prm%nu = lattice_nu(ph) - prm%C66 = lattice_C66(1:6,1:6,ph) + prm%mu = elastic_mu(ph) + prm%nu = elastic_nu(ph) + prm%C66 = elastic_C66(ph) !-------------------------------------------------------------------------------------------------- ! slip related parameters N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray) prm%sum_N_sl = sum(abs(N_sl)) slipActive: if (prm%sum_N_sl > 0) then - prm%P_sl = lattice_SchmidMatrix_slip(N_sl,phase%get_asString('lattice'),& - phase%get_asFloat('c/a',defaultVal=0.0_pReal)) - prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dFloat('h_sl-sl'), & - phase%get_asString('lattice')) - prm%forestProjection = lattice_forestProjection_edge(N_sl,phase%get_asString('lattice'),& - phase%get_asFloat('c/a',defaultVal=0.0_pReal)) + prm%P_sl = lattice_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph)) + prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dFloat('h_sl-sl'),phase_lattice(ph)) + prm%forestProjection = lattice_forestProjection_edge(N_sl,phase_lattice(ph),phase_cOverA(ph)) prm%forestProjection = transpose(prm%forestProjection) - prm%n0_sl = lattice_slip_normal(N_sl,phase%get_asString('lattice'),& - phase%get_asFloat('c/a',defaultVal=0.0_pReal)) + prm%n0_sl = lattice_slip_normal(N_sl,phase_lattice(ph),phase_cOverA(ph)) prm%fccTwinTransNucleation = phase_lattice(ph) == 'cF' .and. (N_sl(1) == 12) - if(prm%fccTwinTransNucleation) prm%fcc_twinNucleationSlipPair = lattice_FCC_TWINNUCLEATIONSLIPPAIR + if (prm%fccTwinTransNucleation) prm%fcc_twinNucleationSlipPair = lattice_FCC_TWINNUCLEATIONSLIPPAIR rho_mob_0 = pl%get_as1dFloat('rho_mob_0', requiredSize=size(N_sl)) rho_dip_0 = pl%get_as1dFloat('rho_dip_0', requiredSize=size(N_sl)) @@ -265,11 +261,9 @@ module function plastic_dislotwin_init() result(myPlasticity) N_tw = pl%get_as1dInt('N_tw', defaultVal=emptyIntArray) prm%sum_N_tw = sum(abs(N_tw)) twinActive: if (prm%sum_N_tw > 0) then - prm%P_tw = lattice_SchmidMatrix_twin(N_tw,phase%get_asString('lattice'),& - phase%get_asFloat('c/a',defaultVal=0.0_pReal)) - prm%h_tw_tw = lattice_interaction_TwinByTwin(N_tw,& - pl%get_as1dFloat('h_tw-tw'), & - phase%get_asString('lattice')) + prm%P_tw = lattice_SchmidMatrix_twin(N_tw,phase_lattice(ph),phase_cOverA(ph)) + prm%h_tw_tw = lattice_interaction_TwinByTwin(N_tw,pl%get_as1dFloat('h_tw-tw'), & + phase_lattice(ph)) prm%b_tw = pl%get_as1dFloat('b_tw', requiredSize=size(N_tw)) prm%t_tw = pl%get_as1dFloat('t_tw', requiredSize=size(N_tw)) @@ -279,11 +273,9 @@ module function plastic_dislotwin_init() result(myPlasticity) prm%L_tw = pl%get_asFloat('L_tw') prm%i_tw = pl%get_asFloat('i_tw') - prm%gamma_char= lattice_characteristicShear_Twin(N_tw,phase%get_asString('lattice'),& - phase%get_asFloat('c/a',defaultVal=0.0_pReal)) + prm%gamma_char= lattice_characteristicShear_Twin(N_tw,phase_lattice(ph),phase_cOverA(ph)) - prm%C66_tw = lattice_C66_twin(N_tw,prm%C66,phase%get_asString('lattice'),& - phase%get_asFloat('c/a',defaultVal=0.0_pReal)) + prm%C66_tw = lattice_C66_twin(N_tw,prm%C66,phase_lattice(ph),phase_cOverA(ph)) if (.not. prm%fccTwinTransNucleation) then prm%dot_N_0_tw = pl%get_as1dFloat('dot_N_0_tw') @@ -324,8 +316,8 @@ module function plastic_dislotwin_init() result(myPlasticity) prm%x_c_tr = pl%get_asFloat('x_c_tr', defaultVal=0.0_pReal) ! ToDo: How to handle that??? prm%L_tr = pl%get_asFloat('L_tr') - prm%h_tr_tr = lattice_interaction_TransByTrans(N_tr,pl%get_as1dFloat('h_tr-tr'), & - phase%get_asString('lattice')) + prm%h_tr_tr = lattice_interaction_TransByTrans(N_tr,pl%get_as1dFloat('h_tr-tr'),& + phase_lattice(ph)) prm%C66_tr = lattice_C66_trans(N_tr,prm%C66,pl%get_asString('lattice_tr'), & 0.0_pReal, & @@ -391,17 +383,15 @@ module function plastic_dislotwin_init() result(myPlasticity) endif slipAndTwinActive: if (prm%sum_N_sl * prm%sum_N_tw > 0) then - prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,N_tw,& - pl%get_as1dFloat('h_sl-tw'), & - phase%get_asString('lattice')) - if (prm%fccTwinTransNucleation .and. size(N_tw) /= 1) extmsg = trim(extmsg)//' interaction_sliptwin' + prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,N_tw,pl%get_as1dFloat('h_sl-tw'), & + phase_lattice(ph)) + if (prm%fccTwinTransNucleation .and. size(N_tw) /= 1) extmsg = trim(extmsg)//' N_tw: nucleation' endif slipAndTwinActive slipAndTransActive: if (prm%sum_N_sl * prm%sum_N_tr > 0) then - prm%h_sl_tr = lattice_interaction_SlipByTrans(N_sl,N_tr,& - pl%get_as1dFloat('h_sl-tr'), & - phase%get_asString('lattice')) - if (prm%fccTwinTransNucleation .and. size(N_tr) /= 1) extmsg = trim(extmsg)//' interaction_sliptrans' + prm%h_sl_tr = lattice_interaction_SlipByTrans(N_sl,N_tr,pl%get_as1dFloat('h_sl-tr'), & + phase_lattice(ph)) + if (prm%fccTwinTransNucleation .and. size(N_tr) /= 1) extmsg = trim(extmsg)//' N_tr: nucleation' endif slipAndTransActive !-------------------------------------------------------------------------------------------------- @@ -491,22 +481,21 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC) real(pReal) :: f_unrotated - associate(prm => param(ph),& - stt => state(ph)) + associate(prm => param(ph), stt => state(ph)) - f_unrotated = 1.0_pReal & - - sum(stt%f_tw(1:prm%sum_N_tw,en)) & - - sum(stt%f_tr(1:prm%sum_N_tr,en)) + f_unrotated = 1.0_pReal & + - sum(stt%f_tw(1:prm%sum_N_tw,en)) & + - sum(stt%f_tr(1:prm%sum_N_tr,en)) - homogenizedC = f_unrotated * prm%C66 - do i=1,prm%sum_N_tw - homogenizedC = homogenizedC & - + stt%f_tw(i,en)*prm%C66_tw(1:6,1:6,i) - enddo - do i=1,prm%sum_N_tr - homogenizedC = homogenizedC & - + stt%f_tr(i,en)*prm%C66_tr(1:6,1:6,i) - enddo + homogenizedC = f_unrotated * prm%C66 + do i=1,prm%sum_N_tw + homogenizedC = homogenizedC & + + stt%f_tw(i,en)*prm%C66_tw(1:6,1:6,i) + enddo + do i=1,prm%sum_N_tr + homogenizedC = homogenizedC & + + stt%f_tr(i,en)*prm%C66_tr(1:6,1:6,i) + enddo end associate @@ -531,7 +520,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en) ddot_gamma_dtau, & tau real(pReal), dimension(param(ph)%sum_N_sl) :: & - dot_gamma_sl,ddot_gamma_dtau_slip + dot_gamma_sl,ddot_gamma_dtau_sl real(pReal), dimension(param(ph)%sum_N_tw) :: & dot_gamma_tw,ddot_gamma_dtau_tw real(pReal), dimension(param(ph)%sum_N_tr) :: & @@ -568,15 +557,15 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en) Lp = 0.0_pReal dLp_dMp = 0.0_pReal - call kinetics_slip(Mp,T,ph,en,dot_gamma_sl,ddot_gamma_dtau_slip) + call kinetics_sl(Mp,T,ph,en,dot_gamma_sl,ddot_gamma_dtau_sl) slipContribution: do i = 1, prm%sum_N_sl Lp = Lp + dot_gamma_sl(i)*prm%P_sl(1:3,1:3,i) forall (k=1:3,l=1:3,m=1:3,n=1:3) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & - + ddot_gamma_dtau_slip(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i) + + ddot_gamma_dtau_sl(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i) enddo slipContribution - call kinetics_twin(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw) + call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw) twinContibution: do i = 1, prm%sum_N_tw Lp = Lp + dot_gamma_tw(i)*prm%P_tw(1:3,1:3,i) forall (k=1:3,l=1:3,m=1:3,n=1:3) & @@ -584,7 +573,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en) + ddot_gamma_dtau_tw(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i) enddo twinContibution - call kinetics_trans(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr,ddot_gamma_dtau_tr) + call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr,ddot_gamma_dtau_tr) transContibution: do i = 1, prm%sum_N_tr Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i) forall (k=1:3,l=1:3,m=1:3,n=1:3) & @@ -609,8 +598,8 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en) StressRatio_p = (abs(tau)/prm%xi_sb)**prm%p_sb dot_gamma_sb = sign(prm%v_sb*exp(-BoltzmannRatio*(1-StressRatio_p)**prm%q_sb), tau) ddot_gamma_dtau = abs(dot_gamma_sb)*BoltzmannRatio* prm%p_sb*prm%q_sb/ prm%xi_sb & - * (abs(tau)/prm%xi_sb)**(prm%p_sb-1.0_pReal) & - * (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal) + * (abs(tau)/prm%xi_sb)**(prm%p_sb-1.0_pReal) & + * (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal) Lp = Lp + dot_gamma_sb * P_sb forall (k=1:3,l=1:3,m=1:3,n=1:3) & @@ -664,7 +653,7 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) - sum(stt%f_tw(1:prm%sum_N_tw,en)) & - sum(stt%f_tr(1:prm%sum_N_tr,en)) - call kinetics_slip(Mp,T,ph,en,dot_gamma_sl) + call kinetics_sl(Mp,T,ph,en,dot_gamma_sl) dot%gamma_sl(:,en) = abs(dot_gamma_sl) rho_dip_distance_min = prm%D_a*prm%b_sl @@ -709,10 +698,10 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) - 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_dip(:,en)*abs(dot_gamma_sl) & - dot_rho_dip_climb - call kinetics_twin(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw) + call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw) dot%f_tw(:,en) = f_unrotated*dot_gamma_tw/prm%gamma_char - call kinetics_trans(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr) + call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr) dot%f_tr(:,en) = f_unrotated*dot_gamma_tr end associate @@ -745,54 +734,52 @@ module subroutine dislotwin_dependentState(T,ph,en) x0 - associate(prm => param(ph),& - stt => state(ph),& - dst => dependentState(ph)) + associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) - sumf_tw = sum(stt%f_tw(1:prm%sum_N_tw,en)) - sumf_tr = sum(stt%f_tr(1:prm%sum_N_tr,en)) + sumf_tw = sum(stt%f_tw(1:prm%sum_N_tw,en)) + sumf_tr = sum(stt%f_tr(1:prm%sum_N_tr,en)) - Gamma = prm%Gamma_sf(1) + prm%Gamma_sf(2) * T + Gamma = prm%Gamma_sf(1) + prm%Gamma_sf(2) * T - !* rescaled volume fraction for topology - f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,en)/prm%t_tw ! this is per system ... - f_over_t_tr = sumf_tr/prm%t_tr ! but this not + !* rescaled volume fraction for topology + f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,en)/prm%t_tw ! this is per system ... + f_over_t_tr = sumf_tr/prm%t_tr ! but this not ! ToDo ...Physically correct, but naming could be adjusted - inv_lambda_sl = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en)))/prm%i_sl - if (prm%sum_N_tw > 0 .and. prm%sum_N_sl > 0) & - inv_lambda_sl = inv_lambda_sl + matmul(prm%h_sl_tw,f_over_t_tw)/(1.0_pReal-sumf_tw) - if (prm%sum_N_tr > 0 .and. prm%sum_N_sl > 0) & - inv_lambda_sl = inv_lambda_sl + matmul(prm%h_sl_tr,f_over_t_tr)/(1.0_pReal-sumf_tr) - dst%Lambda_sl(:,en) = prm%D / (1.0_pReal+prm%D*inv_lambda_sl) + inv_lambda_sl = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en)))/prm%i_sl + if (prm%sum_N_tw > 0 .and. prm%sum_N_sl > 0) & + inv_lambda_sl = inv_lambda_sl + matmul(prm%h_sl_tw,f_over_t_tw)/(1.0_pReal-sumf_tw) + if (prm%sum_N_tr > 0 .and. prm%sum_N_sl > 0) & + inv_lambda_sl = inv_lambda_sl + matmul(prm%h_sl_tr,f_over_t_tr)/(1.0_pReal-sumf_tr) + dst%Lambda_sl(:,en) = prm%D / (1.0_pReal+prm%D*inv_lambda_sl) - inv_lambda_tw_tw = matmul(prm%h_tw_tw,f_over_t_tw)/(1.0_pReal-sumf_tw) - dst%Lambda_tw(:,en) = prm%i_tw*prm%D/(1.0_pReal+prm%D*inv_lambda_tw_tw) + inv_lambda_tw_tw = matmul(prm%h_tw_tw,f_over_t_tw)/(1.0_pReal-sumf_tw) + dst%Lambda_tw(:,en) = prm%i_tw*prm%D/(1.0_pReal+prm%D*inv_lambda_tw_tw) - inv_lambda_tr_tr = matmul(prm%h_tr_tr,f_over_t_tr)/(1.0_pReal-sumf_tr) - dst%Lambda_tr(:,en) = prm%i_tr*prm%D/(1.0_pReal+prm%D*inv_lambda_tr_tr) + inv_lambda_tr_tr = matmul(prm%h_tr_tr,f_over_t_tr)/(1.0_pReal-sumf_tr) + dst%Lambda_tr(:,en) = prm%i_tr*prm%D/(1.0_pReal+prm%D*inv_lambda_tr_tr) - !* threshold stress for dislocation motion - dst%tau_pass(:,en) = prm%mu*prm%b_sl* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en))) + !* threshold stress for dislocation motion + dst%tau_pass(:,en) = prm%mu*prm%b_sl* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en))) - !* threshold stress for growing twin/martensite - if(prm%sum_N_tw == prm%sum_N_sl) & - dst%tau_hat_tw(:,en) = Gamma/(3.0_pReal*prm%b_tw) & - + 3.0_pReal*prm%b_tw*prm%mu/(prm%L_tw*prm%b_sl) ! slip Burgers here correct? - if(prm%sum_N_tr == prm%sum_N_sl) & - dst%tau_hat_tr(:,en) = Gamma/(3.0_pReal*prm%b_tr) & - + 3.0_pReal*prm%b_tr*prm%mu/(prm%L_tr*prm%b_sl) & ! slip Burgers here correct? - + prm%h*prm%delta_G/ (3.0_pReal*prm%b_tr) + !* threshold stress for growing twin/martensite + if(prm%sum_N_tw == prm%sum_N_sl) & + dst%tau_hat_tw(:,en) = Gamma/(3.0_pReal*prm%b_tw) & + + 3.0_pReal*prm%b_tw*prm%mu/(prm%L_tw*prm%b_sl) ! slip Burgers here correct? + if(prm%sum_N_tr == prm%sum_N_sl) & + dst%tau_hat_tr(:,en) = Gamma/(3.0_pReal*prm%b_tr) & + + 3.0_pReal*prm%b_tr*prm%mu/(prm%L_tr*prm%b_sl) & ! slip Burgers here correct? + + prm%h*prm%delta_G/ (3.0_pReal*prm%b_tr) - dst%V_tw(:,en) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,en)**2.0_pReal - dst%V_tr(:,en) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,en)**2.0_pReal + dst%V_tw(:,en) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,en)**2.0_pReal + dst%V_tr(:,en) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,en)**2.0_pReal - x0 = prm%mu*prm%b_tw**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip and is the same for twin and trans - dst%tau_r_tw(:,en) = prm%mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(pi/3.0_pReal)/x0) + x0 = prm%mu*prm%b_tw**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip and is the same for twin and trans + dst%tau_r_tw(:,en) = prm%mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(pi/3.0_pReal)/x0) - x0 = prm%mu*prm%b_tr**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip - dst%tau_r_tr(:,en) = prm%mu*prm%b_tr/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(pi/3.0_pReal)/x0) + x0 = prm%mu*prm%b_tr**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip + dst%tau_r_tr(:,en) = prm%mu*prm%b_tr/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(pi/3.0_pReal)/x0) end associate @@ -857,8 +844,8 @@ end subroutine plastic_dislotwin_results ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end !-------------------------------------------------------------------------------------------------- -pure subroutine kinetics_slip(Mp,T,ph,en, & - dot_gamma_sl,ddot_gamma_dtau_slip,tau_slip) +pure subroutine kinetics_sl(Mp,T,ph,en, & + dot_gamma_sl,ddot_gamma_dtau_sl,tau_sl) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -871,8 +858,8 @@ pure subroutine kinetics_slip(Mp,T,ph,en, & real(pReal), dimension(param(ph)%sum_N_sl), intent(out) :: & dot_gamma_sl real(pReal), dimension(param(ph)%sum_N_sl), optional, intent(out) :: & - ddot_gamma_dtau_slip, & - tau_slip + ddot_gamma_dtau_sl, & + tau_sl real(pReal), dimension(param(ph)%sum_N_sl) :: & ddot_gamma_dtau @@ -891,9 +878,7 @@ pure subroutine kinetics_slip(Mp,T,ph,en, & associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) - do i = 1, prm%sum_N_sl - tau(i) = math_tensordot(Mp,prm%P_sl(1:3,1:3,i)) - enddo + tau = [(math_tensordot(Mp,prm%P_sl(1:3,1:3,i)),i = 1, prm%sum_N_sl)] tau_eff = abs(tau)-dst%tau_pass(:,en) @@ -921,10 +906,10 @@ pure subroutine kinetics_slip(Mp,T,ph,en, & end associate - if(present(ddot_gamma_dtau_slip)) ddot_gamma_dtau_slip = ddot_gamma_dtau - if(present(tau_slip)) tau_slip = tau + if(present(ddot_gamma_dtau_sl)) ddot_gamma_dtau_sl = ddot_gamma_dtau + if(present(tau_sl)) tau_sl = tau -end subroutine kinetics_slip +end subroutine kinetics_sl !-------------------------------------------------------------------------------------------------- @@ -934,8 +919,8 @@ end subroutine kinetics_slip ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- -pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,ph,en,& - dot_gamma_tw,ddot_gamma_dtau_tw) +pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,& + dot_gamma_tw,ddot_gamma_dtau_tw) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -993,7 +978,7 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,ph,en,& if(present(ddot_gamma_dtau_tw)) ddot_gamma_dtau_tw = ddot_gamma_dtau -end subroutine kinetics_twin +end subroutine kinetics_tw !-------------------------------------------------------------------------------------------------- @@ -1003,8 +988,8 @@ end subroutine kinetics_twin ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- -pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,ph,en,& - dot_gamma_tr,ddot_gamma_dtau_tr) +pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,& + dot_gamma_tr,ddot_gamma_dtau_tr) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -1061,6 +1046,6 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,ph,en,& if(present(ddot_gamma_dtau_tr)) ddot_gamma_dtau_tr = ddot_gamma_dtau -end subroutine kinetics_trans +end subroutine kinetics_tr end submodule dislotwin diff --git a/src/phase_mechanical_plastic_isotropic.f90 b/src/phase_mechanical_plastic_isotropic.f90 index 39a619897..8b10ed5ef 100644 --- a/src/phase_mechanical_plastic_isotropic.f90 +++ b/src/phase_mechanical_plastic_isotropic.f90 @@ -170,27 +170,28 @@ module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) integer :: & k, l, m, n + associate(prm => param(ph), stt => state(ph)) - Mp_dev = math_deviatoric33(Mp) - squarenorm_Mp_dev = math_tensordot(Mp_dev,Mp_dev) - norm_Mp_dev = sqrt(squarenorm_Mp_dev) + Mp_dev = math_deviatoric33(Mp) + squarenorm_Mp_dev = math_tensordot(Mp_dev,Mp_dev) + norm_Mp_dev = sqrt(squarenorm_Mp_dev) - if (norm_Mp_dev > 0.0_pReal) then - dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%M*stt%xi(en))) **prm%n + if (norm_Mp_dev > 0.0_pReal) then + dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%M*stt%xi(en))) **prm%n - Lp = dot_gamma * Mp_dev/norm_Mp_dev - forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLp_dMp(k,l,m,n) = (prm%n-1.0_pReal) * Mp_dev(k,l)*Mp_dev(m,n) / squarenorm_Mp_dev - forall (k=1:3,l=1:3) & - dLp_dMp(k,l,k,l) = dLp_dMp(k,l,k,l) + 1.0_pReal - forall (k=1:3,m=1:3) & - dLp_dMp(k,k,m,m) = dLp_dMp(k,k,m,m) - 1.0_pReal/3.0_pReal - dLp_dMp = dot_gamma * dLp_dMp / norm_Mp_dev - else - Lp = 0.0_pReal - dLp_dMp = 0.0_pReal - end if + Lp = dot_gamma * Mp_dev/norm_Mp_dev + forall (k=1:3,l=1:3,m=1:3,n=1:3) & + dLp_dMp(k,l,m,n) = (prm%n-1.0_pReal) * Mp_dev(k,l)*Mp_dev(m,n) / squarenorm_Mp_dev + forall (k=1:3,l=1:3) & + dLp_dMp(k,l,k,l) = dLp_dMp(k,l,k,l) + 1.0_pReal + forall (k=1:3,m=1:3) & + dLp_dMp(k,k,m,m) = dLp_dMp(k,k,m,m) - 1.0_pReal/3.0_pReal + dLp_dMp = dot_gamma * dLp_dMp / norm_Mp_dev + else + Lp = 0.0_pReal + dLp_dMp = 0.0_pReal + end if end associate @@ -218,19 +219,20 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,en) integer :: & k, l, m, n + associate(prm => param(ph), stt => state(ph)) - tr=math_trace33(math_spherical33(Mi)) + tr=math_trace33(math_spherical33(Mi)) - if (prm%dilatation .and. abs(tr) > 0.0_pReal) then ! no stress or J2 plasticity --> Li and its derivative are zero - Li = math_I3 & - * prm%dot_gamma_0 * (3.0_pReal*prm%M*stt%xi(en))**(-prm%n) & - * tr * abs(tr)**(prm%n-1.0_pReal) - forall (k=1:3,l=1:3,m=1:3,n=1:3) dLi_dMi(k,l,m,n) = prm%n / tr * Li(k,l) * math_I3(m,n) - else - Li = 0.0_pReal - dLi_dMi = 0.0_pReal - endif + if (prm%dilatation .and. abs(tr) > 0.0_pReal) then ! no stress or J2 plasticity --> Li and its derivative are zero + Li = math_I3 & + * prm%dot_gamma_0 * (3.0_pReal*prm%M*stt%xi(en))**(-prm%n) & + * tr * abs(tr)**(prm%n-1.0_pReal) + forall (k=1:3,l=1:3,m=1:3,n=1:3) dLi_dMi(k,l,m,n) = prm%n / tr * Li(k,l) * math_I3(m,n) + else + Li = 0.0_pReal + dLi_dMi = 0.0_pReal + endif end associate diff --git a/src/phase_mechanical_plastic_kinehardening.f90 b/src/phase_mechanical_plastic_kinehardening.f90 index 69be33958..880b128a9 100644 --- a/src/phase_mechanical_plastic_kinehardening.f90 +++ b/src/phase_mechanical_plastic_kinehardening.f90 @@ -19,11 +19,11 @@ submodule(phase:plastic) kinehardening xi_inf_f, & xi_inf_b real(pReal), allocatable, dimension(:,:) :: & - interaction_slipslip !< slip resistance from slip activity + h_sl_sl !< slip resistance from slip activity real(pReal), allocatable, dimension(:,:,:) :: & P, & - nonSchmid_pos, & - nonSchmid_neg + P_nS_pos, & + P_nS_neg integer :: & sum_N_sl logical :: & @@ -33,13 +33,14 @@ submodule(phase:plastic) kinehardening end type tParameters type :: tKinehardeningState - real(pReal), pointer, dimension(:,:) :: & !< vectors along NipcMyInstance - crss, & !< critical resolved stress - crss_back, & !< critical resolved back stress - sense, & !< sense of acting shear stress (-1 or +1) - chi0, & !< backstress at last switch of stress sense - gamma0, & !< accumulated shear at last switch of stress sense - accshear !< accumulated (absolute) shear + real(pReal), pointer, dimension(:,:) :: & + xi, & !< resistance against plastic slip + chi, & !< back stress + chi_0, & !< back stress at last switch of stress sense + gamma, & !< accumulated (absolute) shear + gamma_0, & !< accumulated shear at last switch of stress sense + sgn_gamma !< sense of acting shear stress (-1 or +1) + end type tKinehardeningState !-------------------------------------------------------------------------------------------------- @@ -112,21 +113,19 @@ module function plastic_kinehardening_init() result(myPlasticity) N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray) prm%sum_N_sl = sum(abs(N_sl)) slipActive: if (prm%sum_N_sl > 0) then - prm%P = lattice_SchmidMatrix_slip(N_sl,phase%get_asString('lattice'),& - phase%get_asFloat('c/a',defaultVal=0.0_pReal)) + prm%P = lattice_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph)) - if(trim(phase%get_asString('lattice')) == 'cI') then + if (phase_lattice(ph) == 'cI') then a = pl%get_as1dFloat('a_nonSchmid',defaultVal = emptyRealArray) if(size(a) > 0) prm%nonSchmidActive = .true. - prm%nonSchmid_pos = lattice_nonSchmidMatrix(N_sl,a,+1) - prm%nonSchmid_neg = lattice_nonSchmidMatrix(N_sl,a,-1) + prm%P_nS_pos = lattice_nonSchmidMatrix(N_sl,a,+1) + prm%P_nS_neg = lattice_nonSchmidMatrix(N_sl,a,-1) else - prm%nonSchmid_pos = prm%P - prm%nonSchmid_neg = prm%P + prm%P_nS_pos = prm%P + prm%P_nS_neg = prm%P endif - prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(N_sl, & - pl%get_as1dFloat('h_sl-sl'), & - phase%get_asString('lattice')) + prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dFloat('h_sl-sl'), & + phase_lattice(ph)) xi_0 = pl%get_as1dFloat('xi_0', requiredSize=size(N_sl)) prm%xi_inf_f = pl%get_as1dFloat('xi_inf_f', requiredSize=size(N_sl)) @@ -156,18 +155,17 @@ module function plastic_kinehardening_init() result(myPlasticity) if (any(prm%xi_inf_f <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_inf_f' if (any(prm%xi_inf_b <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_inf_b' - !ToDo: Any sensible checks for theta? else slipActive xi_0 = emptyRealArray allocate(prm%xi_inf_f,prm%xi_inf_b,prm%h_0_f,prm%h_inf_f,prm%h_0_b,prm%h_inf_b,source=emptyRealArray) - allocate(prm%interaction_SlipSlip(0,0)) + allocate(prm%h_sl_sl(0,0)) endif slipActive !-------------------------------------------------------------------------------------------------- ! allocate state arrays Nmembers = count(material_phaseID == ph) - sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%sum_N_sl !ToDo: adjust names like in material.yaml - sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl !ToDo: adjust names like in material.yaml + sizeDotState = size(['xi ','chi ', 'gamma']) * prm%sum_N_sl + sizeDeltaState = size(['sgn_gamma', 'chi_0 ', 'gamma_0 ']) * prm%sum_N_sl sizeState = sizeDotState + sizeDeltaState call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState) @@ -176,40 +174,40 @@ module function plastic_kinehardening_init() result(myPlasticity) ! state aliases and initialization startIndex = 1 endIndex = prm%sum_N_sl - stt%crss => plasticState(ph)%state (startIndex:endIndex,:) - stt%crss = spread(xi_0, 2, Nmembers) - dot%crss => plasticState(ph)%dotState(startIndex:endIndex,:) + stt%xi => plasticState(ph)%state (startIndex:endIndex,:) + stt%xi = spread(xi_0, 2, Nmembers) + dot%xi => plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi' startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl - stt%crss_back => plasticState(ph)%state (startIndex:endIndex,:) - dot%crss_back => plasticState(ph)%dotState(startIndex:endIndex,:) + stt%chi => plasticState(ph)%state (startIndex:endIndex,:) + dot%chi => plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl - stt%accshear => plasticState(ph)%state (startIndex:endIndex,:) - dot%accshear => plasticState(ph)%dotState(startIndex:endIndex,:) + stt%gamma => plasticState(ph)%state (startIndex:endIndex,:) + dot%gamma => plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' o = plasticState(ph)%offsetDeltaState startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl - stt%sense => plasticState(ph)%state (startIndex :endIndex ,:) - dlt%sense => plasticState(ph)%deltaState(startIndex-o:endIndex-o,:) + stt%sgn_gamma => plasticState(ph)%state (startIndex :endIndex ,:) + dlt%sgn_gamma => plasticState(ph)%deltaState(startIndex-o:endIndex-o,:) startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl - stt%chi0 => plasticState(ph)%state (startIndex :endIndex ,:) - dlt%chi0 => plasticState(ph)%deltaState(startIndex-o:endIndex-o,:) + stt%chi_0 => plasticState(ph)%state (startIndex :endIndex ,:) + dlt%chi_0 => plasticState(ph)%deltaState(startIndex-o:endIndex-o,:) startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl - stt%gamma0 => plasticState(ph)%state (startIndex :endIndex ,:) - dlt%gamma0 => plasticState(ph)%deltaState(startIndex-o:endIndex-o,:) + stt%gamma_0 => plasticState(ph)%state (startIndex :endIndex ,:) + dlt%gamma_0 => plasticState(ph)%deltaState(startIndex-o:endIndex-o,:) end associate @@ -242,22 +240,22 @@ pure module subroutine kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) integer :: & i,k,l,m,n real(pReal), dimension(param(ph)%sum_N_sl) :: & - gdot_pos,gdot_neg, & - dgdot_dtau_pos,dgdot_dtau_neg + dot_gamma_pos,dot_gamma_neg, & + ddot_gamma_dtau_pos,ddot_gamma_dtau_neg Lp = 0.0_pReal dLp_dMp = 0.0_pReal associate(prm => param(ph)) - call kinetics(Mp,ph,en,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) + call kinetics(Mp,ph,en,dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg) do i = 1, prm%sum_N_sl - Lp = Lp + (gdot_pos(i)+gdot_neg(i))*prm%P(1:3,1:3,i) + Lp = Lp + (dot_gamma_pos(i)+dot_gamma_neg(i))*prm%P(1:3,1:3,i) forall (k=1:3,l=1:3,m=1:3,n=1:3) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & - + dgdot_dtau_pos(i) * prm%P(k,l,i) * prm%nonSchmid_pos(m,n,i) & - + dgdot_dtau_neg(i) * prm%P(k,l,i) * prm%nonSchmid_neg(m,n,i) + + ddot_gamma_dtau_pos(i) * prm%P(k,l,i) * prm%P_nS_pos(m,n,i) & + + ddot_gamma_dtau_neg(i) * prm%P(k,l,i) * prm%P_nS_neg(m,n,i) enddo end associate @@ -279,29 +277,28 @@ module subroutine plastic_kinehardening_dotState(Mp,ph,en) real(pReal) :: & sumGamma real(pReal), dimension(param(ph)%sum_N_sl) :: & - gdot_pos,gdot_neg + dot_gamma_pos,dot_gamma_neg - associate(prm => param(ph), stt => state(ph),& - dot => dotState(ph)) + associate(prm => param(ph), stt => state(ph),dot => dotState(ph)) - call kinetics(Mp,ph,en,gdot_pos,gdot_neg) - dot%accshear(:,en) = abs(gdot_pos+gdot_neg) - sumGamma = sum(stt%accshear(:,en)) + call kinetics(Mp,ph,en,dot_gamma_pos,dot_gamma_neg) + dot%gamma(:,en) = abs(dot_gamma_pos+dot_gamma_neg) + sumGamma = sum(stt%gamma(:,en)) - dot%crss(:,en) = matmul(prm%interaction_SlipSlip,dot%accshear(:,en)) & - * ( prm%h_inf_f & - + (prm%h_0_f - prm%h_inf_f + prm%h_0_f*prm%h_inf_f*sumGamma/prm%xi_inf_f) & - * exp(-sumGamma*prm%h_0_f/prm%xi_inf_f) & - ) + dot%xi(:,en) = matmul(prm%h_sl_sl,dot%gamma(:,en)) & + * ( prm%h_inf_f & + + (prm%h_0_f - prm%h_inf_f + prm%h_0_f*prm%h_inf_f*sumGamma/prm%xi_inf_f) & + * exp(-sumGamma*prm%h_0_f/prm%xi_inf_f) & + ) - dot%crss_back(:,en) = stt%sense(:,en)*dot%accshear(:,en) * & - ( prm%h_inf_b + & - (prm%h_0_b - prm%h_inf_b & - + prm%h_0_b*prm%h_inf_b/(prm%xi_inf_b+stt%chi0(:,en))*(stt%accshear(:,en)-stt%gamma0(:,en))& - ) *exp(-(stt%accshear(:,en)-stt%gamma0(:,en)) *prm%h_0_b/(prm%xi_inf_b+stt%chi0(:,en))) & - ) + dot%chi(:,en) = stt%sgn_gamma(:,en)*dot%gamma(:,en) * & + ( prm%h_inf_b + & + (prm%h_0_b - prm%h_inf_b & + + prm%h_0_b*prm%h_inf_b/(prm%xi_inf_b+stt%chi_0(:,en))*(stt%gamma(:,en)-stt%gamma_0(:,en))& + ) *exp(-(stt%gamma(:,en)-stt%gamma_0(:,en)) *prm%h_0_b/(prm%xi_inf_b+stt%chi_0(:,en))) & + ) end associate @@ -320,28 +317,26 @@ module subroutine plastic_kinehardening_deltaState(Mp,ph,en) en real(pReal), dimension(param(ph)%sum_N_sl) :: & - gdot_pos,gdot_neg, & - sense + dot_gamma_pos,dot_gamma_neg, & + sgn_gamma + associate(prm => param(ph), stt => state(ph), dlt => deltaState(ph)) - call kinetics(Mp,ph,en,gdot_pos,gdot_neg) - sense = merge(state(ph)%sense(:,en), & ! keep existing... - sign(1.0_pReal,gdot_pos+gdot_neg), & ! ...or have a defined - dEq0(gdot_pos+gdot_neg,1e-10_pReal)) ! current sense of shear direction + call kinetics(Mp,ph,en,dot_gamma_pos,dot_gamma_neg) + sgn_gamma = merge(state(ph)%sgn_gamma(:,en), & + sign(1.0_pReal,dot_gamma_pos+dot_gamma_neg), & + dEq0(dot_gamma_pos+dot_gamma_neg,1e-10_pReal)) - -!-------------------------------------------------------------------------------------------------- -! switch in sense of shear? - where(dNeq(sense,stt%sense(:,en),0.1_pReal)) - dlt%sense (:,en) = sense - stt%sense(:,en) ! switch sense - dlt%chi0 (:,en) = abs(stt%crss_back(:,en)) - stt%chi0(:,en) ! remember current backstress magnitude - dlt%gamma0(:,en) = stt%accshear(:,en) - stt%gamma0(:,en) ! remember current accumulated shear - else where - dlt%sense (:,en) = 0.0_pReal - dlt%chi0 (:,en) = 0.0_pReal - dlt%gamma0(:,en) = 0.0_pReal - end where + where(dNeq(sgn_gamma,stt%sgn_gamma(:,en),0.1_pReal)) ! ToDo sgn_gamma*stt%sgn_gamma(:,en)<0 + dlt%sgn_gamma (:,en) = sgn_gamma - stt%sgn_gamma(:,en) + dlt%chi_0 (:,en) = abs(stt%chi(:,en)) - stt%chi_0(:,en) + dlt%gamma_0(:,en) = stt%gamma(:,en) - stt%gamma_0(:,en) + else where + dlt%sgn_gamma (:,en) = 0.0_pReal + dlt%chi_0 (:,en) = 0.0_pReal + dlt%gamma_0(:,en) = 0.0_pReal + end where end associate @@ -362,22 +357,22 @@ module subroutine plastic_kinehardening_results(ph,group) outputsLoop: do o = 1,size(prm%output) select case(trim(prm%output(o))) case ('xi') - if(prm%sum_N_sl>0) call results_writeDataset(stt%crss,group,trim(prm%output(o)), & + if(prm%sum_N_sl>0) call results_writeDataset(stt%xi,group,trim(prm%output(o)), & 'resistance against plastic slip','Pa') - case ('tau_b') - if(prm%sum_N_sl>0) call results_writeDataset(stt%crss_back,group,trim(prm%output(o)), & - 'back stress against plastic slip','Pa') + case ('tau_b') !ToDo: chi + if(prm%sum_N_sl>0) call results_writeDataset(stt%chi,group,trim(prm%output(o)), & + 'back stress','Pa') case ('sgn(gamma)') - if(prm%sum_N_sl>0) call results_writeDataset(stt%sense,group,trim(prm%output(o)), & ! ToDo: could be int + if(prm%sum_N_sl>0) call results_writeDataset(stt%sgn_gamma,group,trim(prm%output(o)), & ! ToDo: could be int 'sense of shear','1') case ('chi_0') - if(prm%sum_N_sl>0) call results_writeDataset(stt%chi0,group,trim(prm%output(o)), & - 'tbd','Pa') + if(prm%sum_N_sl>0) call results_writeDataset(stt%chi_0,group,trim(prm%output(o)), & + 'back stress at last switch of stress sense','Pa') case ('gamma_0') - if(prm%sum_N_sl>0) call results_writeDataset(stt%gamma0,group,trim(prm%output(o)), & - 'tbd','1') + if(prm%sum_N_sl>0) call results_writeDataset(stt%gamma_0,group,trim(prm%output(o)), & + 'plastic shear at last switch of stress sense','1') case ('gamma') - if(prm%sum_N_sl>0) call results_writeDataset(stt%accshear,group,trim(prm%output(o)), & + if(prm%sum_N_sl>0) call results_writeDataset(stt%gamma,group,trim(prm%output(o)), & 'plastic shear','1') end select enddo outputsLoop @@ -394,62 +389,64 @@ end subroutine plastic_kinehardening_results ! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- pure subroutine kinetics(Mp,ph,en, & - gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) + dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress - integer, intent(in) :: & + integer, intent(in) :: & ph, & en real(pReal), intent(out), dimension(param(ph)%sum_N_sl) :: & - gdot_pos, & - gdot_neg - real(pReal), intent(out), optional, dimension(param(ph)%sum_N_sl) :: & - dgdot_dtau_pos, & - dgdot_dtau_neg + dot_gamma_pos, & + dot_gamma_neg + real(pReal), intent(out), dimension(param(ph)%sum_N_sl), optional :: & + ddot_gamma_dtau_pos, & + ddot_gamma_dtau_neg real(pReal), dimension(param(ph)%sum_N_sl) :: & tau_pos, & tau_neg integer :: i + associate(prm => param(ph), stt => state(ph)) - do i = 1, prm%sum_N_sl - tau_pos(i) = math_tensordot(Mp,prm%nonSchmid_pos(1:3,1:3,i)) - stt%crss_back(i,en) - tau_neg(i) = merge(math_tensordot(Mp,prm%nonSchmid_neg(1:3,1:3,i)) - stt%crss_back(i,en), & - 0.0_pReal, prm%nonSchmidActive) - enddo + do i = 1, prm%sum_N_sl + tau_pos(i) = math_tensordot(Mp,prm%P_nS_pos(1:3,1:3,i)) - stt%chi(i,en) + tau_neg(i) = merge(math_tensordot(Mp,prm%P_nS_neg(1:3,1:3,i)) - stt%chi(i,en), & + 0.0_pReal, prm%nonSchmidActive) + enddo - where(dNeq0(tau_pos)) - gdot_pos = prm%dot_gamma_0 * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active - * sign(abs(tau_pos/stt%crss(:,en))**prm%n, tau_pos) - else where - gdot_pos = 0.0_pReal - end where - - where(dNeq0(tau_neg)) - gdot_neg = prm%dot_gamma_0 * 0.5_pReal & ! only used if non-Schmid active, always 1/2 - * sign(abs(tau_neg/stt%crss(:,en))**prm%n, tau_neg) - else where - gdot_neg = 0.0_pReal - end where - - if (present(dgdot_dtau_pos)) then - where(dNeq0(gdot_pos)) - dgdot_dtau_pos = gdot_pos*prm%n/tau_pos + where(dNeq0(tau_pos)) + dot_gamma_pos = prm%dot_gamma_0 * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active + * sign(abs(tau_pos/stt%xi(:,en))**prm%n, tau_pos) else where - dgdot_dtau_pos = 0.0_pReal + dot_gamma_pos = 0.0_pReal end where - endif - if (present(dgdot_dtau_neg)) then - where(dNeq0(gdot_neg)) - dgdot_dtau_neg = gdot_neg*prm%n/tau_neg + + where(dNeq0(tau_neg)) + dot_gamma_neg = prm%dot_gamma_0 * 0.5_pReal & ! only used if non-Schmid active, always 1/2 + * sign(abs(tau_neg/stt%xi(:,en))**prm%n, tau_neg) else where - dgdot_dtau_neg = 0.0_pReal + dot_gamma_neg = 0.0_pReal end where - endif + + if (present(ddot_gamma_dtau_pos)) then + where(dNeq0(dot_gamma_pos)) + ddot_gamma_dtau_pos = dot_gamma_pos*prm%n/tau_pos + else where + ddot_gamma_dtau_pos = 0.0_pReal + end where + endif + if (present(ddot_gamma_dtau_neg)) then + where(dNeq0(dot_gamma_neg)) + ddot_gamma_dtau_neg = dot_gamma_neg*prm%n/tau_neg + else where + ddot_gamma_dtau_neg = 0.0_pReal + end where + endif + end associate end subroutine kinetics diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index fbb092d0e..9510d0b81 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -108,9 +108,9 @@ submodule(phase:plastic) nonlocal forestProjection_Edge, & !< matrix of forest projections of edge dislocations forestProjection_Screw !< matrix of forest projections of screw dislocations real(pReal), dimension(:,:,:), allocatable :: & - Schmid, & !< Schmid contribution - nonSchmid_pos, & - nonSchmid_neg !< combined projection of Schmid and non-Schmid contributions to the resolved shear stress (only for screws) + P_sl, & !< Schmid contribution + P_nS_pos, & + P_nS_neg !< combined projection of Schmid and non-Schmid contributions to the resolved shear stress (only for screws) integer :: & sum_N_sl = 0 integer, dimension(:), allocatable :: & @@ -240,41 +240,35 @@ module function plastic_nonlocal_init() result(myPlasticity) prm%atol_rho = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) - ! This data is read in already in lattice - prm%mu = lattice_mu(ph) - prm%nu = lattice_nu(ph) + prm%mu = elastic_mu(ph) + prm%nu = elastic_nu(ph) ini%N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray) prm%sum_N_sl = sum(abs(ini%N_sl)) slipActive: if (prm%sum_N_sl > 0) then - prm%Schmid = lattice_SchmidMatrix_slip(ini%N_sl,phase%get_asString('lattice'),& - phase%get_asFloat('c/a',defaultVal=0.0_pReal)) + prm%P_sl = lattice_SchmidMatrix_slip(ini%N_sl,phase_lattice(ph), phase_cOverA(ph)) - if(trim(phase%get_asString('lattice')) == 'cI') then + if (phase_lattice(ph) == 'cI') then a = pl%get_as1dFloat('a_nonSchmid',defaultVal = emptyRealArray) if(size(a) > 0) prm%nonSchmidActive = .true. - prm%nonSchmid_pos = lattice_nonSchmidMatrix(ini%N_sl,a,+1) - prm%nonSchmid_neg = lattice_nonSchmidMatrix(ini%N_sl,a,-1) + prm%P_nS_pos = lattice_nonSchmidMatrix(ini%N_sl,a,+1) + prm%P_nS_neg = lattice_nonSchmidMatrix(ini%N_sl,a,-1) else - prm%nonSchmid_pos = prm%Schmid - prm%nonSchmid_neg = prm%Schmid + prm%P_nS_pos = prm%P_sl + prm%P_nS_neg = prm%P_sl endif - prm%h_sl_sl = lattice_interaction_SlipBySlip(ini%N_sl, & - pl%get_as1dFloat('h_sl-sl'), & - phase%get_asString('lattice')) + prm%h_sl_sl = lattice_interaction_SlipBySlip(ini%N_sl,pl%get_as1dFloat('h_sl-sl'), & + phase_lattice(ph)) - prm%forestProjection_edge = lattice_forestProjection_edge (ini%N_sl,phase%get_asString('lattice'),& - phase%get_asFloat('c/a',defaultVal=0.0_pReal)) - prm%forestProjection_screw = lattice_forestProjection_screw(ini%N_sl,phase%get_asString('lattice'),& - phase%get_asFloat('c/a',defaultVal=0.0_pReal)) + prm%forestProjection_edge = lattice_forestProjection_edge (ini%N_sl,phase_lattice(ph),& + phase_cOverA(ph)) + prm%forestProjection_screw = lattice_forestProjection_screw(ini%N_sl,phase_lattice(ph),& + phase_cOverA(ph)) - prm%slip_direction = lattice_slip_direction (ini%N_sl,phase%get_asString('lattice'),& - phase%get_asFloat('c/a',defaultVal=0.0_pReal)) - prm%slip_transverse = lattice_slip_transverse(ini%N_sl,phase%get_asString('lattice'),& - phase%get_asFloat('c/a',defaultVal=0.0_pReal)) - prm%slip_normal = lattice_slip_normal (ini%N_sl,phase%get_asString('lattice'),& - phase%get_asFloat('c/a',defaultVal=0.0_pReal)) + prm%slip_direction = lattice_slip_direction (ini%N_sl,phase_lattice(ph),phase_cOverA(ph)) + prm%slip_transverse = lattice_slip_transverse(ini%N_sl,phase_lattice(ph),phase_cOverA(ph)) + prm%slip_normal = lattice_slip_normal (ini%N_sl,phase_lattice(ph),phase_cOverA(ph)) ! collinear systems (only for octahedral slip systems in fcc) allocate(prm%colinearSystem(prm%sum_N_sl), source = -1) @@ -761,11 +755,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & Mp !< derivative of Lp with respect to Mp integer :: & - ns, & !< short notation for the total number of active slip systems - i, & - j, & - k, & - l, & + i, j, k, l, & t, & !< dislocation type s !< index of my current slip system real(pReal), dimension(param(ph)%sum_N_sl,8) :: & @@ -779,26 +769,24 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & dv_dtauNS !< velocity derivative with respect to the shear stress real(pReal), dimension(param(ph)%sum_N_sl) :: & tau, & !< resolved shear stress including backstress terms - gdotTotal !< shear rate + dot_gamma !< shear rate - associate(prm => param(ph),dst=>microstructure(ph),& - stt=>state(ph)) - ns = prm%sum_N_sl + associate(prm => param(ph),dst=>microstructure(ph),stt=>state(ph)) !*** shortcut to state variables rho = getRho(ph,en) rhoSgl = rho(:,sgl) - do s = 1,ns - tau(s) = math_tensordot(Mp, prm%Schmid(1:3,1:3,s)) + do s = 1,prm%sum_N_sl + tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s)) tauNS(s,1) = tau(s) tauNS(s,2) = tau(s) if (tau(s) > 0.0_pReal) then - tauNS(s,3) = math_tensordot(Mp, +prm%nonSchmid_pos(1:3,1:3,s)) - tauNS(s,4) = math_tensordot(Mp, -prm%nonSchmid_neg(1:3,1:3,s)) + tauNS(s,3) = math_tensordot(Mp, +prm%P_nS_pos(1:3,1:3,s)) + tauNS(s,4) = math_tensordot(Mp, -prm%P_nS_neg(1:3,1:3,s)) else - tauNS(s,3) = math_tensordot(Mp, +prm%nonSchmid_neg(1:3,1:3,s)) - tauNS(s,4) = math_tensordot(Mp, -prm%nonSchmid_pos(1:3,1:3,s)) + tauNS(s,3) = math_tensordot(Mp, +prm%P_nS_neg(1:3,1:3,s)) + tauNS(s,4) = math_tensordot(Mp, -prm%P_nS_pos(1:3,1:3,s)) endif enddo tauNS = tauNS + spread(dst%tau_back(:,en),2,4) @@ -826,22 +814,22 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & stt%v(:,en) = pack(v,.true.) !*** Bauschinger effect - forall (s = 1:ns, t = 5:8, rhoSgl(s,t) * v(s,t-4) < 0.0_pReal) & + forall (s = 1:prm%sum_N_sl, t = 5:8, rhoSgl(s,t) * v(s,t-4) < 0.0_pReal) & rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t)) - gdotTotal = sum(rhoSgl(:,1:4) * v, 2) * prm%b_sl + dot_gamma = sum(rhoSgl(:,1:4) * v, 2) * prm%b_sl Lp = 0.0_pReal dLp_dMp = 0.0_pReal - do s = 1,ns - Lp = Lp + gdotTotal(s) * prm%Schmid(1:3,1:3,s) + do s = 1,prm%sum_N_sl + Lp = Lp + dot_gamma(s) * prm%P_sl(1:3,1:3,s) forall (i=1:3,j=1:3,k=1:3,l=1:3) & dLp_dMp(i,j,k,l) = dLp_dMp(i,j,k,l) & - + prm%Schmid(i,j,s) * prm%Schmid(k,l,s) & + + prm%P_sl(i,j,s) * prm%P_sl(k,l,s) & * sum(rhoSgl(s,1:4) * dv_dtau(s,1:4)) * prm%b_sl(s) & - + prm%Schmid(i,j,s) & - * (+ prm%nonSchmid_pos(k,l,s) * rhoSgl(s,3) * dv_dtauNS(s,3) & - - prm%nonSchmid_neg(k,l,s) * rhoSgl(s,4) * dv_dtauNS(s,4)) * prm%b_sl(s) + + prm%P_sl(i,j,s) & + * (+ prm%P_nS_pos(k,l,s) * rhoSgl(s,3) * dv_dtauNS(s,3) & + - prm%P_nS_neg(k,l,s) * rhoSgl(s,4) * dv_dtauNS(s,4)) * prm%b_sl(s) enddo end associate @@ -861,7 +849,6 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en) en integer :: & - ns, & ! short notation for the total number of active slip systems c, & ! character of dislocation t, & ! type of dislocation s ! index of my current slip system @@ -881,11 +868,10 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en) deltaDUpper ! change in maximum stable dipole distance for edges and screws associate(prm => param(ph),dst => microstructure(ph),del => deltaState(ph)) - ns = prm%sum_N_sl !*** shortcut to state variables - forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en) - forall (s = 1:ns, c = 1:2) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,ph),en) + forall (s = 1:prm%sum_N_sl, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en) + forall (s = 1:prm%sum_N_sl, c = 1:2) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,ph),en) rho = getRho(ph,en) rhoDip = rho(:,dip) @@ -908,7 +894,7 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en) !*** calculate limits for stable dipole height do s = 1,prm%sum_N_sl - tau(s) = math_tensordot(Mp, prm%Schmid(1:3,1:3,s)) +dst%tau_back(s,en) + tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s)) +dst%tau_back(s,en) if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal enddo @@ -926,16 +912,16 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en) !*** dissociation by stress increase deltaRhoDipole2SingleStress = 0.0_pReal - forall (c=1:2, s=1:ns, deltaDUpper(s,c) < 0.0_pReal .and. & + forall (c=1:2, s=1:prm%sum_N_sl, deltaDUpper(s,c) < 0.0_pReal .and. & dNeq0(dUpperOld(s,c) - prm%minDipoleHeight(s,c))) & deltaRhoDipole2SingleStress(s,8+c) = rhoDip(s,c) * deltaDUpper(s,c) & / (dUpperOld(s,c) - prm%minDipoleHeight(s,c)) forall (t=1:4) deltaRhoDipole2SingleStress(:,t) = -0.5_pReal * deltaRhoDipole2SingleStress(:,(t-1)/2+9) - forall (s = 1:ns, c = 1:2) plasticState(ph)%state(iD(s,c,ph),en) = dUpper(s,c) + forall (s = 1:prm%sum_N_sl, c = 1:2) plasticState(ph)%state(iD(s,c,ph),en) = dUpper(s,c) plasticState(ph)%deltaState(:,en) = 0.0_pReal - del%rho(:,en) = reshape(deltaRhoRemobilization + deltaRhoDipole2SingleStress, [10*ns]) + del%rho(:,en) = reshape(deltaRhoRemobilization + deltaRhoDipole2SingleStress, [10*prm%sum_N_sl]) end associate @@ -946,7 +932,7 @@ end subroutine plastic_nonlocal_deltaState !> @brief calculates the rate of change of microstructure !--------------------------------------------------------------------------------------------------- module subroutine nonlocal_dotState(Mp, Temperature,timestep, & - ph,en,ip,el) + ph,en,ip,el) real(pReal), dimension(3,3), intent(in) :: & Mp !< MandelStress @@ -960,7 +946,6 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & el !< current element number integer :: & - ns, & !< short notation for the total number of active slip systems c, & !< character of dislocation t, & !< type of dislocation s !< index of my current slip system @@ -978,7 +963,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & real(pReal), dimension(param(ph)%sum_N_sl,4) :: & v, & !< current dislocation glide velocity v0, & - gdot !< shear rates + dot_gamma !< shear rates real(pReal), dimension(param(ph)%sum_N_sl) :: & tau, & !< current resolved shear stress v_climb !< climb velocity of edge dipoles @@ -994,14 +979,10 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & return endif - associate(prm => param(ph), & - dst => microstructure(ph), & - dot => dotState(ph), & - stt => state(ph)) - ns = prm%sum_N_sl + associate(prm => param(ph), dst => microstructure(ph), dot => dotState(ph), stt => state(ph)) tau = 0.0_pReal - gdot = 0.0_pReal + dot_gamma = 0.0_pReal rho = getRho(ph,en) rhoSgl = rho(:,sgl) @@ -1009,14 +990,14 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & rho0 = getRho0(ph,en) my_rhoSgl0 = rho0(:,sgl) - forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en) - gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4) + forall (s = 1:prm%sum_N_sl, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en) + dot_gamma = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4) ! limits for stable dipole height - do s = 1,ns - tau(s) = math_tensordot(Mp, prm%Schmid(1:3,1:3,s)) + dst%tau_back(s,en) + do s = 1,prm%sum_N_sl + tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s)) + dst%tau_back(s,en) if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal enddo @@ -1035,22 +1016,22 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & ! dislocation multiplication rhoDotMultiplication = 0.0_pReal isBCC: if (phase_lattice(ph) == 'cI') then - forall (s = 1:ns, sum(abs(v(s,1:4))) > 0.0_pReal) - rhoDotMultiplication(s,1:2) = sum(abs(gdot(s,3:4))) / prm%b_sl(s) & ! assuming double-cross-slip of screws to be decisive for multiplication + forall (s = 1:prm%sum_N_sl, sum(abs(v(s,1:4))) > 0.0_pReal) + rhoDotMultiplication(s,1:2) = sum(abs(dot_gamma(s,3:4))) / prm%b_sl(s) & ! assuming double-cross-slip of screws to be decisive for multiplication * sqrt(stt%rho_forest(s,en)) / prm%i_sl(s) ! & ! mean free path ! * 2.0_pReal * sum(abs(v(s,3:4))) / sum(abs(v(s,1:4))) ! ratio of screw to overall velocity determines edge generation - rhoDotMultiplication(s,3:4) = sum(abs(gdot(s,3:4))) /prm%b_sl(s) & ! assuming double-cross-slip of screws to be decisive for multiplication + rhoDotMultiplication(s,3:4) = sum(abs(dot_gamma(s,3:4))) /prm%b_sl(s) & ! assuming double-cross-slip of screws to be decisive for multiplication * sqrt(stt%rho_forest(s,en)) / prm%i_sl(s) ! & ! mean free path ! * 2.0_pReal * sum(abs(v(s,1:2))) / sum(abs(v(s,1:4))) ! ratio of edge to overall velocity determines screw generation endforall else isBCC rhoDotMultiplication(:,1:4) = spread( & - (sum(abs(gdot(:,1:2)),2) * prm%f_ed_mult + sum(abs(gdot(:,3:4)),2)) & + (sum(abs(dot_gamma(:,1:2)),2) * prm%f_ed_mult + sum(abs(dot_gamma(:,3:4)),2)) & * sqrt(stt%rho_forest(:,en)) / prm%i_sl / prm%b_sl, 2, 4) ! eq. 3.26 endif isBCC - forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en) + forall (s = 1:prm%sum_N_sl, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en) !**************************************************************************** @@ -1059,20 +1040,20 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & ! formation by glide do c = 1,2 rhoDotSingle2DipoleGlide(:,2*c-1) = -2.0_pReal * dUpper(:,c) / prm%b_sl & - * ( rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) & ! negative mobile --> positive mobile - + rhoSgl(:,2*c) * abs(gdot(:,2*c-1)) & ! positive mobile --> negative mobile - + abs(rhoSgl(:,2*c+4)) * abs(gdot(:,2*c-1))) ! positive mobile --> negative immobile + * ( rhoSgl(:,2*c-1) * abs(dot_gamma(:,2*c)) & ! negative mobile --> positive mobile + + rhoSgl(:,2*c) * abs(dot_gamma(:,2*c-1)) & ! positive mobile --> negative mobile + + abs(rhoSgl(:,2*c+4)) * abs(dot_gamma(:,2*c-1))) ! positive mobile --> negative immobile rhoDotSingle2DipoleGlide(:,2*c) = -2.0_pReal * dUpper(:,c) / prm%b_sl & - * ( rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) & ! negative mobile --> positive mobile - + rhoSgl(:,2*c) * abs(gdot(:,2*c-1)) & ! positive mobile --> negative mobile - + abs(rhoSgl(:,2*c+3)) * abs(gdot(:,2*c))) ! negative mobile --> positive immobile + * ( rhoSgl(:,2*c-1) * abs(dot_gamma(:,2*c)) & ! negative mobile --> positive mobile + + rhoSgl(:,2*c) * abs(dot_gamma(:,2*c-1)) & ! positive mobile --> negative mobile + + abs(rhoSgl(:,2*c+3)) * abs(dot_gamma(:,2*c))) ! negative mobile --> positive immobile rhoDotSingle2DipoleGlide(:,2*c+3) = -2.0_pReal * dUpper(:,c) / prm%b_sl & - * rhoSgl(:,2*c+3) * abs(gdot(:,2*c)) ! negative mobile --> positive immobile + * rhoSgl(:,2*c+3) * abs(dot_gamma(:,2*c)) ! negative mobile --> positive immobile rhoDotSingle2DipoleGlide(:,2*c+4) = -2.0_pReal * dUpper(:,c) / prm%b_sl & - * rhoSgl(:,2*c+4) * abs(gdot(:,2*c-1)) ! positive mobile --> negative immobile + * rhoSgl(:,2*c+4) * abs(dot_gamma(:,2*c-1)) ! positive mobile --> negative immobile rhoDotSingle2DipoleGlide(:,c+8) = abs(rhoDotSingle2DipoleGlide(:,2*c+3)) & + abs(rhoDotSingle2DipoleGlide(:,2*c+4)) & @@ -1085,13 +1066,13 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & rhoDotAthermalAnnihilation = 0.0_pReal forall (c=1:2) & rhoDotAthermalAnnihilation(:,c+8) = -2.0_pReal * dLower(:,c) / prm%b_sl & - * ( 2.0_pReal * (rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) + rhoSgl(:,2*c) * abs(gdot(:,2*c-1))) & ! was single hitting single - + 2.0_pReal * (abs(rhoSgl(:,2*c+3)) * abs(gdot(:,2*c)) + abs(rhoSgl(:,2*c+4)) * abs(gdot(:,2*c-1))) & ! was single hitting immobile single or was immobile single hit by single - + rhoDip(:,c) * (abs(gdot(:,2*c-1)) + abs(gdot(:,2*c)))) ! single knocks dipole constituent + * ( 2.0_pReal * (rhoSgl(:,2*c-1) * abs(dot_gamma(:,2*c)) + rhoSgl(:,2*c) * abs(dot_gamma(:,2*c-1))) & ! was single hitting single + + 2.0_pReal * (abs(rhoSgl(:,2*c+3)) * abs(dot_gamma(:,2*c)) + abs(rhoSgl(:,2*c+4)) * abs(dot_gamma(:,2*c-1))) & ! was single hitting immobile single or was immobile single hit by single + + rhoDip(:,c) * (abs(dot_gamma(:,2*c-1)) + abs(dot_gamma(:,2*c)))) ! single knocks dipole constituent ! annihilated screw dipoles leave edge jogs behind on the colinear system if (phase_lattice(ph) == 'cF') & - forall (s = 1:ns, prm%colinearSystem(s) > 0) & + forall (s = 1:prm%sum_N_sl, prm%colinearSystem(s) > 0) & rhoDotAthermalAnnihilation(prm%colinearSystem(s),1:2) = - rhoDotAthermalAnnihilation(s,10) & * 0.25_pReal * sqrt(stt%rho_forest(s,en)) * (dUpper(s,2) + dLower(s,2)) * prm%f_ed @@ -1101,7 +1082,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & D_SD = prm%D_0 * exp(-prm%Q_cl / (kB * Temperature)) ! eq. 3.53 v_climb = D_SD * prm%mu * prm%V_at & / (PI * (1.0_pReal-prm%nu) * (dUpper(:,1) + dLower(:,1)) * kB * Temperature) ! eq. 3.54 - forall (s = 1:ns, dUpper(s,1) > dLower(s,1)) & + forall (s = 1:prm%sum_N_sl, dUpper(s,1) > dLower(s,1)) & rhoDotThermalAnnihilation(s,9) = max(- 4.0_pReal * rhoDip(s,1) * v_climb(s) / (dUpper(s,1) - dLower(s,1)), & - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) & - rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have @@ -1124,7 +1105,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & plasticState(ph)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN) else dot%rho(:,en) = pack(rhoDot,.true.) - dot%gamma(:,en) = sum(gdot,2) + dot%gamma(:,en) = sum(dot_gamma,2) endif end associate @@ -1174,7 +1155,7 @@ function rhoDotFlux(timestep,ph,en,ip,el) v, & !< current dislocation glide velocity v0, & neighbor_v0, & !< dislocation glide velocity of enighboring ip - gdot !< shear rates + dot_gamma !< shear rates real(pReal), dimension(3,param(ph)%sum_N_sl,4) :: & m !< direction of dislocation motion real(pReal), dimension(3,3) :: & @@ -1200,7 +1181,7 @@ function rhoDotFlux(timestep,ph,en,ip,el) stt => state(ph)) ns = prm%sum_N_sl - gdot = 0.0_pReal + dot_gamma = 0.0_pReal rho = getRho(ph,en) rhoSgl = rho(:,sgl) @@ -1208,7 +1189,7 @@ function rhoDotFlux(timestep,ph,en,ip,el) my_rhoSgl0 = rho0(:,sgl) forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en) !ToDo: MD: I think we should use state0 here - gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4) + dot_gamma = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4) forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en) @@ -1218,14 +1199,14 @@ function rhoDotFlux(timestep,ph,en,ip,el) if (plasticState(ph)%nonlocal) then !*** check CFL (Courant-Friedrichs-Lewy) condition for flux - if (any( abs(gdot) > 0.0_pReal & ! any active slip system ... + if (any( abs(dot_gamma) > 0.0_pReal & ! any active slip system ... .and. prm%C_CFL * abs(v0) * timestep & > IPvolume(ip,el) / maxval(IParea(:,ip,el)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here) #ifdef DEBUG if (debugConstitutive%extensive) then print'(a,i5,a,i2)', '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip print'(a,e10.3,a,e10.3)', '<< CONST >> velocity is at ', & - maxval(abs(v0), abs(gdot) > 0.0_pReal & + maxval(abs(v0), abs(dot_gamma) > 0.0_pReal & .and. prm%C_CFL * abs(v0) * timestep & > IPvolume(ip,el) / maxval(IParea(:,ip,el))), & ' at a timestep of ',timestep diff --git a/src/phase_mechanical_plastic_phenopowerlaw.f90 b/src/phase_mechanical_plastic_phenopowerlaw.f90 index 45e7a5f14..3bfb360b4 100644 --- a/src/phase_mechanical_plastic_phenopowerlaw.f90 +++ b/src/phase_mechanical_plastic_phenopowerlaw.f90 @@ -33,8 +33,8 @@ submodule(phase:plastic) phenopowerlaw real(pReal), allocatable, dimension(:,:,:) :: & P_sl, & P_tw, & - nonSchmid_pos, & - nonSchmid_neg + P_nS_pos, & + P_nS_neg integer :: & sum_N_sl, & !< total number of active slip system sum_N_tw !< total number of active twin systems @@ -46,10 +46,10 @@ submodule(phase:plastic) phenopowerlaw type :: tPhenopowerlawState real(pReal), pointer, dimension(:,:) :: & - xi_slip, & - xi_twin, & - gamma_slip, & - gamma_twin + xi_sl, & + xi_tw, & + gamma_sl, & + gamma_tw end type tPhenopowerlawState !-------------------------------------------------------------------------------------------------- @@ -115,21 +115,19 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray) prm%sum_N_sl = sum(abs(N_sl)) slipActive: if (prm%sum_N_sl > 0) then - prm%P_sl = lattice_SchmidMatrix_slip(N_sl,phase%get_asString('lattice'),& - phase%get_asFloat('c/a',defaultVal=0.0_pReal)) + prm%P_sl = lattice_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph)) - if(phase%get_asString('lattice') == 'cI') then + if (phase_lattice(ph) == 'cI') then a = pl%get_as1dFloat('a_nonSchmid',defaultVal=emptyRealArray) if(size(a) > 0) prm%nonSchmidActive = .true. - prm%nonSchmid_pos = lattice_nonSchmidMatrix(N_sl,a,+1) - prm%nonSchmid_neg = lattice_nonSchmidMatrix(N_sl,a,-1) + prm%P_nS_pos = lattice_nonSchmidMatrix(N_sl,a,+1) + prm%P_nS_neg = lattice_nonSchmidMatrix(N_sl,a,-1) else - prm%nonSchmid_pos = prm%P_sl - prm%nonSchmid_neg = prm%P_sl + prm%P_nS_pos = prm%P_sl + prm%P_nS_neg = prm%P_sl endif - prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl, & - pl%get_as1dFloat('h_sl-sl'), & - phase%get_asString('lattice')) + prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dFloat('h_sl-sl'), & + phase_lattice(ph)) xi_0_sl = pl%get_as1dFloat('xi_0_sl', requiredSize=size(N_sl)) prm%xi_inf_sl = pl%get_as1dFloat('xi_inf_sl', requiredSize=size(N_sl)) @@ -164,13 +162,11 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) N_tw = pl%get_as1dInt('N_tw', defaultVal=emptyIntArray) prm%sum_N_tw = sum(abs(N_tw)) twinActive: if (prm%sum_N_tw > 0) then - prm%P_tw = lattice_SchmidMatrix_twin(N_tw,phase%get_asString('lattice'),& - phase%get_asFloat('c/a',defaultVal=0.0_pReal)) - prm%h_tw_tw = lattice_interaction_TwinByTwin(N_tw,& - pl%get_as1dFloat('h_tw-tw'), & - phase%get_asString('lattice')) - prm%gamma_char = lattice_characteristicShear_twin(N_tw,phase%get_asString('lattice'),& - phase%get_asFloat('c/a',defaultVal=0.0_pReal)) + prm%P_tw = lattice_SchmidMatrix_twin(N_tw,phase_lattice(ph),phase_cOverA(ph)) + prm%h_tw_tw = lattice_interaction_TwinByTwin(N_tw,pl%get_as1dFloat('h_tw-tw'), & + phase_lattice(ph)) + prm%gamma_char = lattice_characteristicShear_twin(N_tw,phase_lattice(ph),& + phase_cOverA(ph)) xi_0_tw = pl%get_as1dFloat('xi_0_tw',requiredSize=size(N_tw)) @@ -200,12 +196,10 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) ! slip-twin related parameters slipAndTwinActive: if (prm%sum_N_sl > 0 .and. prm%sum_N_tw > 0) then prm%h_0_tw_sl = pl%get_asFloat('h_0_tw-sl') - prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,N_tw,& - pl%get_as1dFloat('h_sl-tw'), & - phase%get_asString('lattice')) - prm%h_tw_sl = lattice_interaction_TwinBySlip(N_tw,N_sl,& - pl%get_as1dFloat('h_tw-sl'), & - phase%get_asString('lattice')) + prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,N_tw,pl%get_as1dFloat('h_sl-tw'), & + phase_lattice(ph)) + prm%h_tw_sl = lattice_interaction_TwinBySlip(N_tw,N_sl,pl%get_as1dFloat('h_tw-sl'), & + phase_lattice(ph)) else slipAndTwinActive allocate(prm%h_sl_tw(prm%sum_N_sl,prm%sum_N_tw)) ! at least one dimension is 0 allocate(prm%h_tw_sl(prm%sum_N_tw,prm%sum_N_sl)) ! at least one dimension is 0 @@ -235,30 +229,30 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) ! state aliases and initialization startIndex = 1 endIndex = prm%sum_N_sl - stt%xi_slip => plasticState(ph)%state (startIndex:endIndex,:) - stt%xi_slip = spread(xi_0_sl, 2, Nmembers) - dot%xi_slip => plasticState(ph)%dotState(startIndex:endIndex,:) + stt%xi_sl => plasticState(ph)%state (startIndex:endIndex,:) + stt%xi_sl = spread(xi_0_sl, 2, Nmembers) + dot%xi_sl => plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi' startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_tw - stt%xi_twin => plasticState(ph)%state (startIndex:endIndex,:) - stt%xi_twin = spread(xi_0_tw, 2, Nmembers) - dot%xi_twin => plasticState(ph)%dotState(startIndex:endIndex,:) + stt%xi_tw => plasticState(ph)%state (startIndex:endIndex,:) + stt%xi_tw = spread(xi_0_tw, 2, Nmembers) + dot%xi_tw => plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl - stt%gamma_slip => plasticState(ph)%state (startIndex:endIndex,:) - dot%gamma_slip => plasticState(ph)%dotState(startIndex:endIndex,:) + stt%gamma_sl => plasticState(ph)%state (startIndex:endIndex,:) + dot%gamma_sl => plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_tw - stt%gamma_twin => plasticState(ph)%state (startIndex:endIndex,:) - dot%gamma_twin => plasticState(ph)%dotState(startIndex:endIndex,:) + stt%gamma_tw => plasticState(ph)%state (startIndex:endIndex,:) + dot%gamma_tw => plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) end associate @@ -293,31 +287,31 @@ pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) integer :: & i,k,l,m,n real(pReal), dimension(param(ph)%sum_N_sl) :: & - gdot_slip_pos,gdot_slip_neg, & - dgdot_dtauslip_pos,dgdot_dtauslip_neg + dot_gamma_sl_pos,dot_gamma_sl_neg, & + ddot_gamma_dtau_sl_pos,ddot_gamma_dtau_sl_neg real(pReal), dimension(param(ph)%sum_N_tw) :: & - gdot_twin,dgdot_dtautwin + dot_gamma_tw,ddot_gamma_dtau_tw Lp = 0.0_pReal dLp_dMp = 0.0_pReal associate(prm => param(ph)) - call kinetics_slip(Mp,ph,en,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg) + call kinetics_sl(Mp,ph,en,dot_gamma_sl_pos,dot_gamma_sl_neg,ddot_gamma_dtau_sl_pos,ddot_gamma_dtau_sl_neg) slipSystems: do i = 1, prm%sum_N_sl - Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%P_sl(1:3,1:3,i) + Lp = Lp + (dot_gamma_sl_pos(i)+dot_gamma_sl_neg(i))*prm%P_sl(1:3,1:3,i) forall (k=1:3,l=1:3,m=1:3,n=1:3) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & - + dgdot_dtauslip_pos(i) * prm%P_sl(k,l,i) * prm%nonSchmid_pos(m,n,i) & - + dgdot_dtauslip_neg(i) * prm%P_sl(k,l,i) * prm%nonSchmid_neg(m,n,i) + + ddot_gamma_dtau_sl_pos(i) * prm%P_sl(k,l,i) * prm%P_nS_pos(m,n,i) & + + ddot_gamma_dtau_sl_neg(i) * prm%P_sl(k,l,i) * prm%P_nS_neg(m,n,i) enddo slipSystems - call kinetics_twin(Mp,ph,en,gdot_twin,dgdot_dtautwin) + call kinetics_tw(Mp,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw) twinSystems: do i = 1, prm%sum_N_tw - Lp = Lp + gdot_twin(i)*prm%P_tw(1:3,1:3,i) + Lp = Lp + dot_gamma_tw(i)*prm%P_tw(1:3,1:3,i) forall (k=1:3,l=1:3,m=1:3,n=1:3) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & - + dgdot_dtautwin(i)*prm%P_tw(k,l,i)*prm%P_tw(m,n,i) + + ddot_gamma_dtau_tw(i)*prm%P_tw(k,l,i)*prm%P_tw(m,n,i) enddo twinSystems end associate @@ -337,46 +331,32 @@ module subroutine phenopowerlaw_dotState(Mp,ph,en) en real(pReal) :: & - c_SlipSlip,c_TwinSlip,c_TwinTwin, & - xi_slip_sat_offset,& - sumGamma,sumF + xi_sl_sat_offset,& + sumF real(pReal), dimension(param(ph)%sum_N_sl) :: & - left_SlipSlip,right_SlipSlip, & - gdot_slip_pos,gdot_slip_neg + dot_gamma_sl_pos,dot_gamma_sl_neg, & + right_SlipSlip - associate(prm => param(ph), stt => state(ph), & - dot => dotState(ph)) - sumGamma = sum(stt%gamma_slip(:,en)) - sumF = sum(stt%gamma_twin(:,en)/prm%gamma_char) + associate(prm => param(ph), stt => state(ph), dot => dotState(ph)) -!-------------------------------------------------------------------------------------------------- -! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices - c_SlipSlip = prm%h_0_sl_sl * (1.0_pReal + prm%c_1*sumF** prm%c_2) - c_TwinSlip = prm%h_0_tw_sl * sumGamma**prm%c_3 - c_TwinTwin = prm%h_0_tw_tw * sumF**prm%c_4 + call kinetics_sl(Mp,ph,en,dot_gamma_sl_pos,dot_gamma_sl_neg) + dot%gamma_sl(:,en) = abs(dot_gamma_sl_pos+dot_gamma_sl_neg) + call kinetics_tw(Mp,ph,en,dot%gamma_tw(:,en)) -!-------------------------------------------------------------------------------------------------- -! calculate left and right vectors - left_SlipSlip = 1.0_pReal + prm%h_int - xi_slip_sat_offset = prm%f_sat_sl_tw*sqrt(sumF) - right_SlipSlip = sign(abs(1.0_pReal-stt%xi_slip(:,en) / (prm%xi_inf_sl+xi_slip_sat_offset)) **prm%a_sl, & - 1.0_pReal-stt%xi_slip(:,en) / (prm%xi_inf_sl+xi_slip_sat_offset)) + sumF = sum(stt%gamma_tw(:,en)/prm%gamma_char) + xi_sl_sat_offset = prm%f_sat_sl_tw*sqrt(sumF) + right_SlipSlip = sign(abs(1.0_pReal-stt%xi_sl(:,en) / (prm%xi_inf_sl+xi_sl_sat_offset))**prm%a_sl, & + 1.0_pReal-stt%xi_sl(:,en) / (prm%xi_inf_sl+xi_sl_sat_offset)) -!-------------------------------------------------------------------------------------------------- -! shear rates - call kinetics_slip(Mp,ph,en,gdot_slip_pos,gdot_slip_neg) - dot%gamma_slip(:,en) = abs(gdot_slip_pos+gdot_slip_neg) - call kinetics_twin(Mp,ph,en,dot%gamma_twin(:,en)) + dot%xi_sl(:,en) = prm%h_0_sl_sl * (1.0_pReal + prm%c_1*sumF** prm%c_2) * (1.0_pReal + prm%h_int) & + * matmul(prm%h_sl_sl,dot%gamma_sl(:,en)*right_SlipSlip) & + + matmul(prm%h_sl_tw,dot%gamma_tw(:,en)) -!-------------------------------------------------------------------------------------------------- -! hardening - dot%xi_slip(:,en) = c_SlipSlip * left_SlipSlip * & - matmul(prm%h_sl_sl,dot%gamma_slip(:,en)*right_SlipSlip) & - + matmul(prm%h_sl_tw,dot%gamma_twin(:,en)) + dot%xi_tw(:,en) = prm%h_0_tw_sl * sum(stt%gamma_sl(:,en))**prm%c_3 & + * matmul(prm%h_tw_sl,dot%gamma_sl(:,en)) & + + prm%h_0_tw_tw * sumF**prm%c_4 * matmul(prm%h_tw_tw,dot%gamma_tw(:,en)) - dot%xi_twin(:,en) = c_TwinSlip * matmul(prm%h_tw_sl,dot%gamma_slip(:,en)) & - + c_TwinTwin * matmul(prm%h_tw_tw,dot%gamma_twin(:,en)) end associate end subroutine phenopowerlaw_dotState @@ -397,17 +377,17 @@ module subroutine plastic_phenopowerlaw_results(ph,group) select case(trim(prm%output(o))) case('xi_sl') - if(prm%sum_N_sl>0) call results_writeDataset(stt%xi_slip,group,trim(prm%output(o)), & + if(prm%sum_N_sl>0) call results_writeDataset(stt%xi_sl,group,trim(prm%output(o)), & 'resistance against plastic slip','Pa') case('gamma_sl') - if(prm%sum_N_sl>0) call results_writeDataset(stt%gamma_slip,group,trim(prm%output(o)), & + if(prm%sum_N_sl>0) call results_writeDataset(stt%gamma_sl,group,trim(prm%output(o)), & 'plastic shear','1') case('xi_tw') - if(prm%sum_N_tw>0) call results_writeDataset(stt%xi_twin,group,trim(prm%output(o)), & + if(prm%sum_N_tw>0) call results_writeDataset(stt%xi_tw,group,trim(prm%output(o)), & 'resistance against twinning','Pa') case('gamma_tw') - if(prm%sum_N_tw>0) call results_writeDataset(stt%gamma_twin,group,trim(prm%output(o)), & + if(prm%sum_N_tw>0) call results_writeDataset(stt%gamma_tw,group,trim(prm%output(o)), & 'twinning shear','1') end select @@ -424,8 +404,8 @@ end subroutine plastic_phenopowerlaw_results ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- -pure subroutine kinetics_slip(Mp,ph,en, & - gdot_slip_pos,gdot_slip_neg,dgdot_dtau_slip_pos,dgdot_dtau_slip_neg) +pure subroutine kinetics_sl(Mp,ph,en, & + dot_gamma_sl_pos,dot_gamma_sl_neg,ddot_gamma_dtau_sl_pos,ddot_gamma_dtau_sl_neg) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -434,56 +414,57 @@ pure subroutine kinetics_slip(Mp,ph,en, & en real(pReal), intent(out), dimension(param(ph)%sum_N_sl) :: & - gdot_slip_pos, & - gdot_slip_neg + dot_gamma_sl_pos, & + dot_gamma_sl_neg real(pReal), intent(out), optional, dimension(param(ph)%sum_N_sl) :: & - dgdot_dtau_slip_pos, & - dgdot_dtau_slip_neg + ddot_gamma_dtau_sl_pos, & + ddot_gamma_dtau_sl_neg real(pReal), dimension(param(ph)%sum_N_sl) :: & - tau_slip_pos, & - tau_slip_neg + tau_sl_pos, & + tau_sl_neg integer :: i associate(prm => param(ph), stt => state(ph)) - do i = 1, prm%sum_N_sl - tau_slip_pos(i) = math_tensordot(Mp,prm%nonSchmid_pos(1:3,1:3,i)) - tau_slip_neg(i) = merge(math_tensordot(Mp,prm%nonSchmid_neg(1:3,1:3,i)), & + do i = 1, prm%sum_N_sl + tau_sl_pos(i) = math_tensordot(Mp,prm%P_nS_pos(1:3,1:3,i)) + tau_sl_neg(i) = merge(math_tensordot(Mp,prm%P_nS_neg(1:3,1:3,i)), & 0.0_pReal, prm%nonSchmidActive) - enddo + enddo - where(dNeq0(tau_slip_pos)) - gdot_slip_pos = prm%dot_gamma_0_sl * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active - * sign(abs(tau_slip_pos/stt%xi_slip(:,en))**prm%n_sl, tau_slip_pos) - else where - gdot_slip_pos = 0.0_pReal - end where - - where(dNeq0(tau_slip_neg)) - gdot_slip_neg = prm%dot_gamma_0_sl * 0.5_pReal & ! only used if non-Schmid active, always 1/2 - * sign(abs(tau_slip_neg/stt%xi_slip(:,en))**prm%n_sl, tau_slip_neg) - else where - gdot_slip_neg = 0.0_pReal - end where - - if (present(dgdot_dtau_slip_pos)) then - where(dNeq0(gdot_slip_pos)) - dgdot_dtau_slip_pos = gdot_slip_pos*prm%n_sl/tau_slip_pos + where(dNeq0(tau_sl_pos)) + dot_gamma_sl_pos = prm%dot_gamma_0_sl * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active + * sign(abs(tau_sl_pos/stt%xi_sl(:,en))**prm%n_sl, tau_sl_pos) else where - dgdot_dtau_slip_pos = 0.0_pReal + dot_gamma_sl_pos = 0.0_pReal end where - endif - if (present(dgdot_dtau_slip_neg)) then - where(dNeq0(gdot_slip_neg)) - dgdot_dtau_slip_neg = gdot_slip_neg*prm%n_sl/tau_slip_neg + + where(dNeq0(tau_sl_neg)) + dot_gamma_sl_neg = prm%dot_gamma_0_sl * 0.5_pReal & ! only used if non-Schmid active, always 1/2 + * sign(abs(tau_sl_neg/stt%xi_sl(:,en))**prm%n_sl, tau_sl_neg) else where - dgdot_dtau_slip_neg = 0.0_pReal + dot_gamma_sl_neg = 0.0_pReal end where - endif + + if (present(ddot_gamma_dtau_sl_pos)) then + where(dNeq0(dot_gamma_sl_pos)) + ddot_gamma_dtau_sl_pos = dot_gamma_sl_pos*prm%n_sl/tau_sl_pos + else where + ddot_gamma_dtau_sl_pos = 0.0_pReal + end where + endif + if (present(ddot_gamma_dtau_sl_neg)) then + where(dNeq0(dot_gamma_sl_neg)) + ddot_gamma_dtau_sl_neg = dot_gamma_sl_neg*prm%n_sl/tau_sl_neg + else where + ddot_gamma_dtau_sl_neg = 0.0_pReal + end where + endif + end associate -end subroutine kinetics_slip +end subroutine kinetics_sl !-------------------------------------------------------------------------------------------------- @@ -493,8 +474,8 @@ end subroutine kinetics_slip ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- -pure subroutine kinetics_twin(Mp,ph,en,& - gdot_twin,dgdot_dtau_twin) +pure subroutine kinetics_tw(Mp,ph,en,& + dot_gamma_tw,ddot_gamma_dtau_tw) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -503,37 +484,36 @@ pure subroutine kinetics_twin(Mp,ph,en,& en real(pReal), dimension(param(ph)%sum_N_tw), intent(out) :: & - gdot_twin + dot_gamma_tw real(pReal), dimension(param(ph)%sum_N_tw), intent(out), optional :: & - dgdot_dtau_twin + ddot_gamma_dtau_tw real(pReal), dimension(param(ph)%sum_N_tw) :: & - tau_twin + tau_tw integer :: i + associate(prm => param(ph), stt => state(ph)) - do i = 1, prm%sum_N_tw - tau_twin(i) = math_tensordot(Mp,prm%P_tw(1:3,1:3,i)) - enddo + tau_tw = [(math_tensordot(Mp,prm%P_tw(1:3,1:3,i)),i=1,prm%sum_N_tw)] - where(tau_twin > 0.0_pReal) - gdot_twin = (1.0_pReal-sum(stt%gamma_twin(:,en)/prm%gamma_char)) & ! only twin in untwinned volume fraction - * prm%dot_gamma_0_tw*(abs(tau_twin)/stt%xi_twin(:,en))**prm%n_tw - else where - gdot_twin = 0.0_pReal - end where - - if (present(dgdot_dtau_twin)) then - where(dNeq0(gdot_twin)) - dgdot_dtau_twin = gdot_twin*prm%n_tw/tau_twin + where(tau_tw > 0.0_pReal) + dot_gamma_tw = (1.0_pReal-sum(stt%gamma_tw(:,en)/prm%gamma_char)) & ! only twin in untwinned volume fraction + * prm%dot_gamma_0_tw*(abs(tau_tw)/stt%xi_tw(:,en))**prm%n_tw else where - dgdot_dtau_twin = 0.0_pReal + dot_gamma_tw = 0.0_pReal end where - endif + + if (present(ddot_gamma_dtau_tw)) then + where(dNeq0(dot_gamma_tw)) + ddot_gamma_dtau_tw = dot_gamma_tw*prm%n_tw/tau_tw + else where + ddot_gamma_dtau_tw = 0.0_pReal + end where + endif end associate -end subroutine kinetics_twin +end subroutine kinetics_tw end submodule phenopowerlaw diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90 index 8d2915ef7..9a75263ca 100644 --- a/src/phase_thermal.f90 +++ b/src/phase_thermal.f90 @@ -103,7 +103,7 @@ module subroutine thermal_init(phases) param(ph)%C_p = thermal%get_asFloat('C_p',defaultVal=0.0_pReal) ! ToDo: make mandatory? param(ph)%K(1,1) = thermal%get_asFloat('K_11',defaultVal=0.0_pReal) ! ToDo: make mandatory? param(ph)%K(3,3) = thermal%get_asFloat('K_33',defaultVal=0.0_pReal) ! ToDo: depends on symmtery - param(ph)%K = lattice_applyLatticeSymmetry33(param(ph)%K,phase_lattice(ph)) + param(ph)%K = lattice_symmetrize_33(param(ph)%K,phase_lattice(ph)) sources => thermal%get('source',defaultVal=emptyList) thermal_Nsources(ph) = sources%length @@ -216,7 +216,7 @@ module function phase_K_T(co,ce) result(K) end function phase_K_T -module function thermal_stress(Delta_t,ph,en) result(converged_) ! ?? why is this called "stress" when it seems closer to "updateState" ?? +module function phase_thermal_constitutive(Delta_t,ph,en) result(converged_) real(pReal), intent(in) :: Delta_t integer, intent(in) :: ph, en @@ -225,7 +225,7 @@ module function thermal_stress(Delta_t,ph,en) result(converged_) ! ?? converged_ = .not. integrateThermalState(Delta_t,ph,en) -end function thermal_stress +end function phase_thermal_constitutive !-------------------------------------------------------------------------------------------------- diff --git a/src/prec.f90 b/src/prec.f90 index 7613cfe46..620647b82 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -38,16 +38,14 @@ module prec sizeDotState = 0, & !< size of dot state, i.e. state(1:sizeDot) follows time evolution by dotState rates offsetDeltaState = 0, & !< index offset of delta state sizeDeltaState = 0 !< size of delta state, i.e. state(offset+1:offset+sizeDelta) follows time evolution by deltaState increments - ! http://stackoverflow.com/questions/3948210 - real(pReal), pointer, dimension(:), contiguous :: & + real(pReal), allocatable, dimension(:) :: & atol - real(pReal), pointer, dimension(:,:), contiguous :: & ! a pointer is needed here because we might point to state/doState. However, they will never point to something, but are rather allocated and, hence, contiguous + ! http://stackoverflow.com/questions/3948210 + real(pReal), pointer, dimension(:,:), contiguous :: & !< is basically an allocatable+target, but in a type needs to be pointer state0, & state, & !< state dotState, & !< rate of state change deltaState !< increment of state change - real(pReal), allocatable, dimension(:,:) :: & - subState0 end type type, extends(tState) :: tPlasticState @@ -61,12 +59,9 @@ module prec real(pReal), private, parameter :: PREAL_EPSILON = epsilon(0.0_pReal) !< minimum positive number such that 1.0 + EPSILON /= 1.0. real(pReal), private, parameter :: PREAL_MIN = tiny(0.0_pReal) !< smallest normalized floating point number - integer, dimension(0), parameter :: & - emptyIntArray = [integer::] - real(pReal), dimension(0), parameter :: & - emptyRealArray = [real(pReal)::] - character(len=pStringLen), dimension(0), parameter :: & - emptyStringArray = [character(len=pStringLen)::] + integer, dimension(0), parameter :: emptyIntArray = [integer::] + real(pReal), dimension(0), parameter :: emptyRealArray = [real(pReal)::] + character(len=pStringLen), dimension(0), parameter :: emptyStringArray = [character(len=pStringLen)::] private :: & selfTest