From 8995faae8ba0c0a4277ca1d2dc67cb018ca998e9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 3 May 2020 10:49:36 +0200 Subject: [PATCH 1/9] WIP: part of dotState should be in deltaState --- src/constitutive_plastic_nonlocal.f90 | 404 ++++++++++++++++++++++++++ 1 file changed, 404 insertions(+) diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index 77f1556f6..153e12dd1 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -1344,6 +1344,410 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, & end subroutine plastic_nonlocal_dotState +!--------------------------------------------------------------------------------------------------- +!> @brief calculates the rate of change of microstructure +!--------------------------------------------------------------------------------------------------- +subroutine plastic_nonlocal_dotState2(Mp, F, Fp, Temperature,timestep, & + instance,of,ip,el) + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< MandelStress + real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: & + F, & !< elastic deformation gradient + Fp !< plastic deformation gradient + real(pReal), intent(in) :: & + Temperature, & !< temperature + timestep !< substepped crystallite time increment + integer, intent(in) :: & + instance, & + of, & + ip, & !< current integration point + el !< current element number + + integer :: & + ph, & + neighbor_instance, & !< instance of my neighbor's plasticity + ns, & !< short notation for the total number of active slip systems + c, & !< character of dislocation + n, & !< index of my current neighbor + neighbor_el, & !< element number of my neighbor + neighbor_ip, & !< integration point of my neighbor + neighbor_n, & !< neighbor index pointing to me when looking from my neighbor + opposite_neighbor, & !< index of my opposite neighbor + opposite_ip, & !< ip of my opposite neighbor + opposite_el, & !< element index of my opposite neighbor + opposite_n, & !< neighbor index pointing to me when looking from my opposite neighbor + t, & !< type of dislocation + no,& !< neighbor offset shortcut + np,& !< neighbor phase shortcut + topp, & !< type of dislocation with opposite sign to t + s !< index of my current slip system + real(pReal), dimension(param(instance)%sum_N_sl,10) :: & + rho, & + rho0, & !< dislocation density at beginning of time step + rhoDot, & !< density evolution + rhoDotMultiplication, & !< density evolution by multiplication + rhoDotFlux, & !< density evolution by flux + rhoDotSingle2DipoleGlide, & !< density evolution by dipole formation (by glide) + rhoDotAthermalAnnihilation, & !< density evolution by athermal annihilation + rhoDotThermalAnnihilation !< density evolution by thermal annihilation + real(pReal), dimension(param(instance)%sum_N_sl,8) :: & + rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles) + neighbor_rhoSgl0, & !< current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles) + my_rhoSgl0 !< single dislocation densities of central ip (positive/negative screw and edge without dipoles) + real(pReal), dimension(param(instance)%sum_N_sl,4) :: & + v, & !< current dislocation glide velocity + v0, & + neighbor_v0, & !< dislocation glide velocity of enighboring ip + gdot !< shear rates + real(pReal), dimension(param(instance)%sum_N_sl) :: & + tau, & !< current resolved shear stress + vClimb !< climb velocity of edge dipoles + real(pReal), dimension(param(instance)%sum_N_sl,2) :: & + rhoDip, & !< current dipole dislocation densities (screw and edge dipoles) + dLower, & !< minimum stable dipole distance for edges and screws + dUpper !< current maximum stable dipole distance for edges and screws + real(pReal), dimension(3,param(instance)%sum_N_sl,4) :: & + m !< direction of dislocation motion + real(pReal), dimension(3,3) :: & + my_F, & !< my total deformation gradient + neighbor_F, & !< total deformation gradient of my neighbor + my_Fe, & !< my elastic deformation gradient + neighbor_Fe, & !< elastic deformation gradient of my neighbor + Favg !< average total deformation gradient of me and my neighbor + real(pReal), dimension(3) :: & + normal_neighbor2me, & !< interface normal pointing from my neighbor to me in neighbor's lattice configuration + normal_neighbor2me_defConf, & !< interface normal pointing from my neighbor to me in shared deformed configuration + normal_me2neighbor, & !< interface normal pointing from me to my neighbor in my lattice configuration + normal_me2neighbor_defConf !< interface normal pointing from me to my neighbor in shared deformed configuration + real(pReal) :: & + area, & !< area of the current interface + transmissivity, & !< overall transmissivity of dislocation flux to neighboring material point + lineLength, & !< dislocation line length leaving the current interface + selfDiffusion !< self diffusion + + ph = material_phaseAt(1,el) + if (timestep <= 0.0_pReal) then + plasticState(ph)%dotState = 0.0_pReal + return + endif + + associate(prm => param(instance), & + dst => microstructure(instance), & + dot => dotState(instance), & + stt => state(instance)) + ns = prm%sum_N_sl + + tau = 0.0_pReal + gdot = 0.0_pReal + + rho = getRho(instance,of,ip,el) + rhoSgl = rho(:,sgl) + rhoDip = rho(:,dip) + rho0 = getRho0(instance,of,ip,el) + my_rhoSgl0 = rho0(:,sgl) + + forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,instance),of) + gdot = rhoSgl(:,1:4) * v * spread(prm%burgers,2,4) + +#ifdef DEBUG + if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0 & + .and. ((debug_e == el .and. debug_i == ip)& + .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0 )) then + write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> rho / 1/m^2', rhoSgl, rhoDip + write(6,'(a,/,4(12x,12(e12.5,1x),/))') '<< CONST >> gdot / 1/s',gdot + endif +#endif + + !**************************************************************************** + !*** 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,of) + if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal + enddo + + dLower = prm%minDipoleHeight + dUpper(:,1) = prm%mu * prm%burgers/(8.0_pReal * PI * (1.0_pReal - prm%nu) * abs(tau)) + dUpper(:,2) = prm%mu * prm%burgers/(4.0_pReal * PI * abs(tau)) + + where(dNeq0(sqrt(sum(abs(rho(:,edg)),2)))) & + dUpper(:,1) = min(1.0_pReal/sqrt(sum(abs(rho(:,edg)),2)),dUpper(:,1)) + where(dNeq0(sqrt(sum(abs(rho(:,scr)),2)))) & + dUpper(:,2) = min(1.0_pReal/sqrt(sum(abs(rho(:,scr)),2)),dUpper(:,2)) + + dUpper = max(dUpper,dLower) + + !**************************************************************************** + !*** dislocation multiplication + rhoDotMultiplication = 0.0_pReal + isBCC: if (lattice_structure(ph) == LATTICE_bcc_ID) 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%burgers(s) & ! assuming double-cross-slip of screws to be decisive for multiplication + * sqrt(stt%rho_forest(s,of)) / prm%lambda0(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%burgers(s) & ! assuming double-cross-slip of screws to be decisive for multiplication + * sqrt(stt%rho_forest(s,of)) / prm%lambda0(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%fEdgeMultiplication + sum(abs(gdot(:,3:4)),2)) & + * sqrt(stt%rho_forest(:,of)) / prm%lambda0 / prm%burgers, 2, 4) + endif isBCC + + forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,instance),of) + + !**************************************************************************** + !*** calculate dislocation fluxes (only for nonlocal plasticity) + rhoDotFlux = 0.0_pReal + if (.not. phase_localPlasticity(material_phaseAt(1,el))) then + + !*** check CFL (Courant-Friedrichs-Lewy) condition for flux + if (any( abs(gdot) > 0.0_pReal & ! any active slip system ... + .and. prm%CFLfactor * 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 (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0) then + write(6,'(a,i5,a,i2)') '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip + write(6,'(a,e10.3,a,e10.3)') '<< CONST >> velocity is at ', & + maxval(abs(v0), abs(gdot) > 0.0_pReal & + .and. prm%CFLfactor * abs(v0) * timestep & + > IPvolume(ip,el) / maxval(IParea(:,ip,el))), & + ' at a timestep of ',timestep + write(6,'(a)') '<< CONST >> enforcing cutback !!!' + endif +#endif + plasticState(ph)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN) ! -> return NaN and, hence, enforce cutback + return + endif + + + !*** be aware of the definition of slip_transverse = slip_direction x slip_normal !!! + !*** opposite sign to our t vector in the (s,t,n) triplet !!! + + m(1:3,:,1) = prm%slip_direction + m(1:3,:,2) = -prm%slip_direction + m(1:3,:,3) = -prm%slip_transverse + m(1:3,:,4) = prm%slip_transverse + + my_F = F(1:3,1:3,1,ip,el) + my_Fe = matmul(my_F, math_inv33(Fp(1:3,1:3,1,ip,el))) + + neighbors: do n = 1,nIPneighbors + + neighbor_el = IPneighborhood(1,n,ip,el) + neighbor_ip = IPneighborhood(2,n,ip,el) + neighbor_n = IPneighborhood(3,n,ip,el) + np = material_phaseAt(1,neighbor_el) + no = material_phasememberAt(1,neighbor_ip,neighbor_el) + + opposite_neighbor = n + mod(n,2) - mod(n+1,2) + opposite_el = IPneighborhood(1,opposite_neighbor,ip,el) + opposite_ip = IPneighborhood(2,opposite_neighbor,ip,el) + opposite_n = IPneighborhood(3,opposite_neighbor,ip,el) + + if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient + neighbor_instance = phase_plasticityInstance(material_phaseAt(1,neighbor_el)) + neighbor_F = F(1:3,1:3,1,neighbor_ip,neighbor_el) + neighbor_Fe = matmul(neighbor_F, math_inv33(Fp(1:3,1:3,1,neighbor_ip,neighbor_el))) + Favg = 0.5_pReal * (my_F + neighbor_F) + else ! if no neighbor, take my value as average + Favg = my_F + endif + + neighbor_v0 = 0.0_pReal ! needed for check of sign change in flux density below + + !* FLUX FROM MY NEIGHBOR TO ME + !* This is only considered, if I have a neighbor of nonlocal plasticity + !* (also nonlocal constitutive law with local properties) that is at least a little bit + !* compatible. + !* If it's not at all compatible, no flux is arriving, because everything is dammed in front of + !* my neighbor's interface. + !* The entering flux from my neighbor will be distributed on my slip systems according to the + !* compatibility + if (neighbor_n > 0) then + if (phase_plasticity(material_phaseAt(1,neighbor_el)) == PLASTICITY_NONLOCAL_ID .and. & + any(compatibility(:,:,:,n,ip,el) > 0.0_pReal)) then + + forall (s = 1:ns, t = 1:4) + neighbor_v0(s,t) = plasticState(np)%state0(iV (s,t,neighbor_instance),no) + neighbor_rhoSgl0(s,t) = max(plasticState(np)%state0(iRhoU(s,t,neighbor_instance),no),0.0_pReal) + endforall + + where (neighbor_rhoSgl0 * IPvolume(neighbor_ip,neighbor_el) ** 0.667_pReal < prm%significantN & + .or. neighbor_rhoSgl0 < prm%significantRho) & + neighbor_rhoSgl0 = 0.0_pReal + normal_neighbor2me_defConf = math_det33(Favg) * matmul(math_inv33(transpose(Favg)), & + IPareaNormal(1:3,neighbor_n,neighbor_ip,neighbor_el)) ! normal of the interface in (average) deformed configuration (pointing neighbor => me) + normal_neighbor2me = matmul(transpose(neighbor_Fe), normal_neighbor2me_defConf) & + / math_det33(neighbor_Fe) ! interface normal in the lattice configuration of my neighbor + area = IParea(neighbor_n,neighbor_ip,neighbor_el) * norm2(normal_neighbor2me) + normal_neighbor2me = normal_neighbor2me / norm2(normal_neighbor2me) ! normalize the surface normal to unit length + do s = 1,ns + do t = 1,4 + c = (t + 1) / 2 + topp = t + mod(t,2) - mod(t+1,2) + if (neighbor_v0(s,t) * math_inner(m(1:3,s,t), normal_neighbor2me) > 0.0_pReal & ! flux from my neighbor to me == entering flux for me + .and. v0(s,t) * neighbor_v0(s,t) >= 0.0_pReal ) then ! ... only if no sign change in flux density + lineLength = neighbor_rhoSgl0(s,t) * neighbor_v0(s,t) & + * math_inner(m(1:3,s,t), normal_neighbor2me) * area ! positive line length that wants to enter through this interface + where (compatibility(c,:,s,n,ip,el) > 0.0_pReal) & + rhoDotFlux(:,t) = rhoDotFlux(1:ns,t) & + + lineLength/IPvolume(ip,el)*compatibility(c,:,s,n,ip,el)**2.0_pReal ! transferring to equally signed mobile dislocation type + where (compatibility(c,:,s,n,ip,el) < 0.0_pReal) & + rhoDotFlux(:,topp) = rhoDotFlux(:,topp) & + + lineLength/IPvolume(ip,el)*compatibility(c,:,s,n,ip,el)**2.0_pReal ! transferring to opposite signed mobile dislocation type + + endif + enddo + enddo + endif; endif + + + !* FLUX FROM ME TO MY NEIGHBOR + !* This is not considered, if my opposite neighbor has a different constitutive law than nonlocal (still considered for nonlocal law with local properties). + !* Then, we assume, that the opposite(!) neighbor sends an equal amount of dislocations to me. + !* So the net flux in the direction of my neighbor is equal to zero: + !* leaving flux to neighbor == entering flux from opposite neighbor + !* In case of reduced transmissivity, part of the leaving flux is stored as dead dislocation density. + !* That means for an interface of zero transmissivity the leaving flux is fully converted to dead dislocations. + if (opposite_n > 0) then + if (phase_plasticity(material_phaseAt(1,opposite_el)) == PLASTICITY_NONLOCAL_ID) then + + normal_me2neighbor_defConf = math_det33(Favg) & + * matmul(math_inv33(transpose(Favg)),IPareaNormal(1:3,n,ip,el)) ! normal of the interface in (average) deformed configuration (pointing me => neighbor) + normal_me2neighbor = matmul(transpose(my_Fe), normal_me2neighbor_defConf) & + / math_det33(my_Fe) ! interface normal in my lattice configuration + area = IParea(n,ip,el) * norm2(normal_me2neighbor) + normal_me2neighbor = normal_me2neighbor / norm2(normal_me2neighbor) ! normalize the surface normal to unit length + do s = 1,ns + do t = 1,4 + c = (t + 1) / 2 + if (v0(s,t) * math_inner(m(1:3,s,t), normal_me2neighbor) > 0.0_pReal ) then ! flux from me to my neighbor == leaving flux for me (might also be a pure flux from my mobile density to dead density if interface not at all transmissive) + if (v0(s,t) * neighbor_v0(s,t) >= 0.0_pReal) then ! no sign change in flux density + transmissivity = sum(compatibility(c,:,s,n,ip,el)**2.0_pReal) ! overall transmissivity from this slip system to my neighbor + else ! sign change in flux density means sign change in stress which does not allow for dislocations to arive at the neighbor + transmissivity = 0.0_pReal + endif + lineLength = my_rhoSgl0(s,t) * v0(s,t) & + * math_inner(m(1:3,s,t), normal_me2neighbor) * area ! positive line length of mobiles that wants to leave through this interface + rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / IPvolume(ip,el) ! subtract dislocation flux from current type + rhoDotFlux(s,t+4) = rhoDotFlux(s,t+4) & + + lineLength / IPvolume(ip,el) * (1.0_pReal - transmissivity) & + * sign(1.0_pReal, v0(s,t)) ! dislocation flux that is not able to leave through interface (because of low transmissivity) will remain as immobile single density at the material point + endif + enddo + enddo + endif; endif + + enddo neighbors + endif + + + + !**************************************************************************** + !*** calculate dipole formation and annihilation + + !*** formation by glide + do c = 1,2 + rhoDotSingle2DipoleGlide(:,2*c-1) = -2.0_pReal * dUpper(:,c) / prm%burgers & + * ( 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 + + rhoDotSingle2DipoleGlide(:,2*c) = -2.0_pReal * dUpper(:,c) / prm%burgers & + * ( 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 + + rhoDotSingle2DipoleGlide(:,2*c+3) = -2.0_pReal * dUpper(:,c) / prm%burgers & + * rhoSgl(:,2*c+3) * abs(gdot(:,2*c)) ! negative mobile --> positive immobile + + rhoDotSingle2DipoleGlide(:,2*c+4) = -2.0_pReal * dUpper(:,c) / prm%burgers & + * rhoSgl(:,2*c+4) * abs(gdot(:,2*c-1)) ! positive mobile --> negative immobile + + rhoDotSingle2DipoleGlide(:,c+8) = abs(rhoDotSingle2DipoleGlide(:,2*c+3)) & + + abs(rhoDotSingle2DipoleGlide(:,2*c+4)) & + - rhoDotSingle2DipoleGlide(:,2*c-1) & + - rhoDotSingle2DipoleGlide(:,2*c) + enddo + + + !*** athermal annihilation + rhoDotAthermalAnnihilation = 0.0_pReal + forall (c=1:2) & + rhoDotAthermalAnnihilation(:,c+8) = -2.0_pReal * dLower(:,c) / prm%burgers & + * ( 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 + + ! annihilated screw dipoles leave edge jogs behind on the colinear system + if (lattice_structure(ph) == LATTICE_fcc_ID) & + forall (s = 1:ns, prm%colinearSystem(s) > 0) & + rhoDotAthermalAnnihilation(prm%colinearSystem(s),1:2) = - rhoDotAthermalAnnihilation(s,10) & + * 0.25_pReal * sqrt(stt%rho_forest(s,of)) * (dUpper(s,2) + dLower(s,2)) * prm%edgeJogFactor + + + !*** thermally activated annihilation of edge dipoles by climb + rhoDotThermalAnnihilation = 0.0_pReal + selfDiffusion = prm%Dsd0 * exp(-prm%selfDiffusionEnergy / (kB * Temperature)) + vClimb = prm%atomicVolume * selfDiffusion * prm%mu & + / ( kB * Temperature * PI * (1.0_pReal-prm%nu) * (dUpper(:,1) + dLower(:,1))) + forall (s = 1:ns, dUpper(s,1) > dLower(s,1)) & + rhoDotThermalAnnihilation(s,9) = max(- 4.0_pReal * rhoDip(s,1) * vClimb(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 + + rhoDot = rhoDotFlux & + + rhoDotMultiplication & + + rhoDotSingle2DipoleGlide & + + rhoDotAthermalAnnihilation & + + rhoDotThermalAnnihilation + +#ifdef DEBUG + if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0 & + .and. ((debug_e == el .and. debug_i == ip)& + .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0 )) then + write(6,'(a,/,4(12x,12(e12.5,1x),/))') '<< CONST >> dislocation multiplication', & + rhoDotMultiplication(:,1:4) * timestep + write(6,'(a,/,8(12x,12(e12.5,1x),/))') '<< CONST >> dislocation flux', & + rhoDotFlux(:,1:8) * timestep + write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> dipole formation by glide', & + rhoDotSingle2DipoleGlide * timestep + write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> athermal dipole annihilation', & + rhoDotAthermalAnnihilation * timestep + write(6,'(a,/,2(12x,12(e12.5,1x),/))') '<< CONST >> thermally activated dipole annihilation', & + rhoDotThermalAnnihilation(:,9:10) * timestep + write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> total density change', & + rhoDot * timestep + write(6,'(a,/,10(12x,12(f12.5,1x),/))') '<< CONST >> relative density change', & + rhoDot(:,1:8) * timestep / (abs(stt%rho(:,sgl))+1.0e-10), & + rhoDot(:,9:10) * timestep / (stt%rho(:,dip)+1.0e-10) + write(6,*) + endif +#endif + + + if ( any(rho(:,mob) + rhoDot(:,1:4) * timestep < -prm%atol_rho) & + .or. any(rho(:,dip) + rhoDot(:,9:10) * timestep < -prm%atol_rho)) then +#ifdef DEBUG + if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0) then + write(6,'(a,i5,a,i2)') '<< CONST >> evolution rate leads to negative density at el ',el,' ip ',ip + write(6,'(a)') '<< CONST >> enforcing cutback !!!' + endif +#endif + plasticState(ph)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN) + else + dot%rho(:,of) = pack(rhoDot,.true.) + dot%gamma(:,of) = sum(gdot,2) + endif + + end associate + +end subroutine plastic_nonlocal_dotState2 + + !-------------------------------------------------------------------------------------------------- !> @brief Compatibility update !> @detail Compatibility is defined as normalized product of signed cosine of the angle between the slip From 6f627892bc9871b73d0ce66c19bdb4cb1365cad7 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 3 May 2020 10:54:57 +0200 Subject: [PATCH 2/9] not needed --- src/constitutive_plastic_nonlocal.f90 | 170 +------------------------- 1 file changed, 4 insertions(+), 166 deletions(-) diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index 153e12dd1..e80f2ceab 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -1347,16 +1347,13 @@ end subroutine plastic_nonlocal_dotState !--------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !--------------------------------------------------------------------------------------------------- -subroutine plastic_nonlocal_dotState2(Mp, F, Fp, Temperature,timestep, & - instance,of,ip,el) +subroutine plastic_nonlocal_dotState2(F, Fp,timestep, & + instance,of,ip,el) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< MandelStress real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: & F, & !< elastic deformation gradient Fp !< plastic deformation gradient real(pReal), intent(in) :: & - Temperature, & !< temperature timestep !< substepped crystallite time increment integer, intent(in) :: & instance, & @@ -1385,12 +1382,7 @@ subroutine plastic_nonlocal_dotState2(Mp, F, Fp, Temperature,timestep, & real(pReal), dimension(param(instance)%sum_N_sl,10) :: & rho, & rho0, & !< dislocation density at beginning of time step - rhoDot, & !< density evolution - rhoDotMultiplication, & !< density evolution by multiplication - rhoDotFlux, & !< density evolution by flux - rhoDotSingle2DipoleGlide, & !< density evolution by dipole formation (by glide) - rhoDotAthermalAnnihilation, & !< density evolution by athermal annihilation - rhoDotThermalAnnihilation !< density evolution by thermal annihilation + rhoDotFlux !< density evolution by flux real(pReal), dimension(param(instance)%sum_N_sl,8) :: & rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles) neighbor_rhoSgl0, & !< current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles) @@ -1400,13 +1392,6 @@ subroutine plastic_nonlocal_dotState2(Mp, F, Fp, Temperature,timestep, & v0, & neighbor_v0, & !< dislocation glide velocity of enighboring ip gdot !< shear rates - real(pReal), dimension(param(instance)%sum_N_sl) :: & - tau, & !< current resolved shear stress - vClimb !< climb velocity of edge dipoles - real(pReal), dimension(param(instance)%sum_N_sl,2) :: & - rhoDip, & !< current dipole dislocation densities (screw and edge dipoles) - dLower, & !< minimum stable dipole distance for edges and screws - dUpper !< current maximum stable dipole distance for edges and screws real(pReal), dimension(3,param(instance)%sum_N_sl,4) :: & m !< direction of dislocation motion real(pReal), dimension(3,3) :: & @@ -1423,8 +1408,7 @@ subroutine plastic_nonlocal_dotState2(Mp, F, Fp, Temperature,timestep, & real(pReal) :: & area, & !< area of the current interface transmissivity, & !< overall transmissivity of dislocation flux to neighboring material point - lineLength, & !< dislocation line length leaving the current interface - selfDiffusion !< self diffusion + lineLength !< dislocation line length leaving the current interface ph = material_phaseAt(1,el) if (timestep <= 0.0_pReal) then @@ -1438,63 +1422,16 @@ subroutine plastic_nonlocal_dotState2(Mp, F, Fp, Temperature,timestep, & stt => state(instance)) ns = prm%sum_N_sl - tau = 0.0_pReal gdot = 0.0_pReal rho = getRho(instance,of,ip,el) rhoSgl = rho(:,sgl) - rhoDip = rho(:,dip) rho0 = getRho0(instance,of,ip,el) my_rhoSgl0 = rho0(:,sgl) forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,instance),of) gdot = rhoSgl(:,1:4) * v * spread(prm%burgers,2,4) -#ifdef DEBUG - if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0 & - .and. ((debug_e == el .and. debug_i == ip)& - .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0 )) then - write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> rho / 1/m^2', rhoSgl, rhoDip - write(6,'(a,/,4(12x,12(e12.5,1x),/))') '<< CONST >> gdot / 1/s',gdot - endif -#endif - - !**************************************************************************** - !*** 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,of) - if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal - enddo - - dLower = prm%minDipoleHeight - dUpper(:,1) = prm%mu * prm%burgers/(8.0_pReal * PI * (1.0_pReal - prm%nu) * abs(tau)) - dUpper(:,2) = prm%mu * prm%burgers/(4.0_pReal * PI * abs(tau)) - - where(dNeq0(sqrt(sum(abs(rho(:,edg)),2)))) & - dUpper(:,1) = min(1.0_pReal/sqrt(sum(abs(rho(:,edg)),2)),dUpper(:,1)) - where(dNeq0(sqrt(sum(abs(rho(:,scr)),2)))) & - dUpper(:,2) = min(1.0_pReal/sqrt(sum(abs(rho(:,scr)),2)),dUpper(:,2)) - - dUpper = max(dUpper,dLower) - - !**************************************************************************** - !*** dislocation multiplication - rhoDotMultiplication = 0.0_pReal - isBCC: if (lattice_structure(ph) == LATTICE_bcc_ID) 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%burgers(s) & ! assuming double-cross-slip of screws to be decisive for multiplication - * sqrt(stt%rho_forest(s,of)) / prm%lambda0(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%burgers(s) & ! assuming double-cross-slip of screws to be decisive for multiplication - * sqrt(stt%rho_forest(s,of)) / prm%lambda0(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%fEdgeMultiplication + sum(abs(gdot(:,3:4)),2)) & - * sqrt(stt%rho_forest(:,of)) / prm%lambda0 / prm%burgers, 2, 4) - endif isBCC forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,instance),of) @@ -1644,105 +1581,6 @@ subroutine plastic_nonlocal_dotState2(Mp, F, Fp, Temperature,timestep, & enddo neighbors endif - - - !**************************************************************************** - !*** calculate dipole formation and annihilation - - !*** formation by glide - do c = 1,2 - rhoDotSingle2DipoleGlide(:,2*c-1) = -2.0_pReal * dUpper(:,c) / prm%burgers & - * ( 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 - - rhoDotSingle2DipoleGlide(:,2*c) = -2.0_pReal * dUpper(:,c) / prm%burgers & - * ( 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 - - rhoDotSingle2DipoleGlide(:,2*c+3) = -2.0_pReal * dUpper(:,c) / prm%burgers & - * rhoSgl(:,2*c+3) * abs(gdot(:,2*c)) ! negative mobile --> positive immobile - - rhoDotSingle2DipoleGlide(:,2*c+4) = -2.0_pReal * dUpper(:,c) / prm%burgers & - * rhoSgl(:,2*c+4) * abs(gdot(:,2*c-1)) ! positive mobile --> negative immobile - - rhoDotSingle2DipoleGlide(:,c+8) = abs(rhoDotSingle2DipoleGlide(:,2*c+3)) & - + abs(rhoDotSingle2DipoleGlide(:,2*c+4)) & - - rhoDotSingle2DipoleGlide(:,2*c-1) & - - rhoDotSingle2DipoleGlide(:,2*c) - enddo - - - !*** athermal annihilation - rhoDotAthermalAnnihilation = 0.0_pReal - forall (c=1:2) & - rhoDotAthermalAnnihilation(:,c+8) = -2.0_pReal * dLower(:,c) / prm%burgers & - * ( 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 - - ! annihilated screw dipoles leave edge jogs behind on the colinear system - if (lattice_structure(ph) == LATTICE_fcc_ID) & - forall (s = 1:ns, prm%colinearSystem(s) > 0) & - rhoDotAthermalAnnihilation(prm%colinearSystem(s),1:2) = - rhoDotAthermalAnnihilation(s,10) & - * 0.25_pReal * sqrt(stt%rho_forest(s,of)) * (dUpper(s,2) + dLower(s,2)) * prm%edgeJogFactor - - - !*** thermally activated annihilation of edge dipoles by climb - rhoDotThermalAnnihilation = 0.0_pReal - selfDiffusion = prm%Dsd0 * exp(-prm%selfDiffusionEnergy / (kB * Temperature)) - vClimb = prm%atomicVolume * selfDiffusion * prm%mu & - / ( kB * Temperature * PI * (1.0_pReal-prm%nu) * (dUpper(:,1) + dLower(:,1))) - forall (s = 1:ns, dUpper(s,1) > dLower(s,1)) & - rhoDotThermalAnnihilation(s,9) = max(- 4.0_pReal * rhoDip(s,1) * vClimb(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 - - rhoDot = rhoDotFlux & - + rhoDotMultiplication & - + rhoDotSingle2DipoleGlide & - + rhoDotAthermalAnnihilation & - + rhoDotThermalAnnihilation - -#ifdef DEBUG - if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0 & - .and. ((debug_e == el .and. debug_i == ip)& - .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0 )) then - write(6,'(a,/,4(12x,12(e12.5,1x),/))') '<< CONST >> dislocation multiplication', & - rhoDotMultiplication(:,1:4) * timestep - write(6,'(a,/,8(12x,12(e12.5,1x),/))') '<< CONST >> dislocation flux', & - rhoDotFlux(:,1:8) * timestep - write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> dipole formation by glide', & - rhoDotSingle2DipoleGlide * timestep - write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> athermal dipole annihilation', & - rhoDotAthermalAnnihilation * timestep - write(6,'(a,/,2(12x,12(e12.5,1x),/))') '<< CONST >> thermally activated dipole annihilation', & - rhoDotThermalAnnihilation(:,9:10) * timestep - write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> total density change', & - rhoDot * timestep - write(6,'(a,/,10(12x,12(f12.5,1x),/))') '<< CONST >> relative density change', & - rhoDot(:,1:8) * timestep / (abs(stt%rho(:,sgl))+1.0e-10), & - rhoDot(:,9:10) * timestep / (stt%rho(:,dip)+1.0e-10) - write(6,*) - endif -#endif - - - if ( any(rho(:,mob) + rhoDot(:,1:4) * timestep < -prm%atol_rho) & - .or. any(rho(:,dip) + rhoDot(:,9:10) * timestep < -prm%atol_rho)) then -#ifdef DEBUG - if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0) then - write(6,'(a,i5,a,i2)') '<< CONST >> evolution rate leads to negative density at el ',el,' ip ',ip - write(6,'(a)') '<< CONST >> enforcing cutback !!!' - endif -#endif - plasticState(ph)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN) - else - dot%rho(:,of) = pack(rhoDot,.true.) - dot%gamma(:,of) = sum(gdot,2) - endif - end associate end subroutine plastic_nonlocal_dotState2 From 53de95798a38e3b85f90b4168a38c74b2671b320 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 3 May 2020 11:49:16 +0200 Subject: [PATCH 3/9] separated local and nonlocal part of dotState --- src/constitutive_plastic_nonlocal.f90 | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index e80f2ceab..928800e1d 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -1295,7 +1295,7 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, & - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) & - rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have - rhoDot = rhoDotFlux & + rhoDot = plastic_nonlocal_dotState2(F,Fp,timestep, instance,of,ip,el) & + rhoDotMultiplication & + rhoDotSingle2DipoleGlide & + rhoDotAthermalAnnihilation & @@ -1347,8 +1347,7 @@ end subroutine plastic_nonlocal_dotState !--------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !--------------------------------------------------------------------------------------------------- -subroutine plastic_nonlocal_dotState2(F, Fp,timestep, & - instance,of,ip,el) +function plastic_nonlocal_dotState2(F,Fp,timestep, instance,of,ip,el) result(rhoDotFlux) real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: & F, & !< elastic deformation gradient @@ -1583,7 +1582,7 @@ subroutine plastic_nonlocal_dotState2(F, Fp,timestep, & end associate -end subroutine plastic_nonlocal_dotState2 +end function plastic_nonlocal_dotState2 !-------------------------------------------------------------------------------------------------- From 2a1badb548ccd5979461ab17344ceeec0a4a7976 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 3 May 2020 12:04:57 +0200 Subject: [PATCH 4/9] not needed --- src/constitutive_plastic_nonlocal.f90 | 180 +------------------------- 1 file changed, 1 insertion(+), 179 deletions(-) diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index 928800e1d..08d14a5ce 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -962,39 +962,24 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, & integer :: & ph, & - neighbor_instance, & !< instance of my neighbor's plasticity ns, & !< short notation for the total number of active slip systems c, & !< character of dislocation - n, & !< index of my current neighbor - neighbor_el, & !< element number of my neighbor - neighbor_ip, & !< integration point of my neighbor - neighbor_n, & !< neighbor index pointing to me when looking from my neighbor - opposite_neighbor, & !< index of my opposite neighbor - opposite_ip, & !< ip of my opposite neighbor - opposite_el, & !< element index of my opposite neighbor - opposite_n, & !< neighbor index pointing to me when looking from my opposite neighbor t, & !< type of dislocation - no,& !< neighbor offset shortcut - np,& !< neighbor phase shortcut - topp, & !< type of dislocation with opposite sign to t s !< index of my current slip system real(pReal), dimension(param(instance)%sum_N_sl,10) :: & rho, & rho0, & !< dislocation density at beginning of time step rhoDot, & !< density evolution rhoDotMultiplication, & !< density evolution by multiplication - rhoDotFlux, & !< density evolution by flux rhoDotSingle2DipoleGlide, & !< density evolution by dipole formation (by glide) rhoDotAthermalAnnihilation, & !< density evolution by athermal annihilation rhoDotThermalAnnihilation !< density evolution by thermal annihilation real(pReal), dimension(param(instance)%sum_N_sl,8) :: & rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles) - neighbor_rhoSgl0, & !< current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles) my_rhoSgl0 !< single dislocation densities of central ip (positive/negative screw and edge without dipoles) real(pReal), dimension(param(instance)%sum_N_sl,4) :: & v, & !< current dislocation glide velocity v0, & - neighbor_v0, & !< dislocation glide velocity of enighboring ip gdot !< shear rates real(pReal), dimension(param(instance)%sum_N_sl) :: & tau, & !< current resolved shear stress @@ -1003,23 +988,7 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, & rhoDip, & !< current dipole dislocation densities (screw and edge dipoles) dLower, & !< minimum stable dipole distance for edges and screws dUpper !< current maximum stable dipole distance for edges and screws - real(pReal), dimension(3,param(instance)%sum_N_sl,4) :: & - m !< direction of dislocation motion - real(pReal), dimension(3,3) :: & - my_F, & !< my total deformation gradient - neighbor_F, & !< total deformation gradient of my neighbor - my_Fe, & !< my elastic deformation gradient - neighbor_Fe, & !< elastic deformation gradient of my neighbor - Favg !< average total deformation gradient of me and my neighbor - real(pReal), dimension(3) :: & - normal_neighbor2me, & !< interface normal pointing from my neighbor to me in neighbor's lattice configuration - normal_neighbor2me_defConf, & !< interface normal pointing from my neighbor to me in shared deformed configuration - normal_me2neighbor, & !< interface normal pointing from me to my neighbor in my lattice configuration - normal_me2neighbor_defConf !< interface normal pointing from me to my neighbor in shared deformed configuration real(pReal) :: & - area, & !< area of the current interface - transmissivity, & !< overall transmissivity of dislocation flux to neighboring material point - lineLength, & !< dislocation line length leaving the current interface selfDiffusion !< self diffusion ph = material_phaseAt(1,el) @@ -1094,153 +1063,6 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, & forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,instance),of) - !**************************************************************************** - !*** calculate dislocation fluxes (only for nonlocal plasticity) - rhoDotFlux = 0.0_pReal - if (.not. phase_localPlasticity(material_phaseAt(1,el))) then - - !*** check CFL (Courant-Friedrichs-Lewy) condition for flux - if (any( abs(gdot) > 0.0_pReal & ! any active slip system ... - .and. prm%CFLfactor * 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 (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0) then - write(6,'(a,i5,a,i2)') '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip - write(6,'(a,e10.3,a,e10.3)') '<< CONST >> velocity is at ', & - maxval(abs(v0), abs(gdot) > 0.0_pReal & - .and. prm%CFLfactor * abs(v0) * timestep & - > IPvolume(ip,el) / maxval(IParea(:,ip,el))), & - ' at a timestep of ',timestep - write(6,'(a)') '<< CONST >> enforcing cutback !!!' - endif -#endif - plasticState(ph)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN) ! -> return NaN and, hence, enforce cutback - return - endif - - - !*** be aware of the definition of slip_transverse = slip_direction x slip_normal !!! - !*** opposite sign to our t vector in the (s,t,n) triplet !!! - - m(1:3,:,1) = prm%slip_direction - m(1:3,:,2) = -prm%slip_direction - m(1:3,:,3) = -prm%slip_transverse - m(1:3,:,4) = prm%slip_transverse - - my_F = F(1:3,1:3,1,ip,el) - my_Fe = matmul(my_F, math_inv33(Fp(1:3,1:3,1,ip,el))) - - neighbors: do n = 1,nIPneighbors - - neighbor_el = IPneighborhood(1,n,ip,el) - neighbor_ip = IPneighborhood(2,n,ip,el) - neighbor_n = IPneighborhood(3,n,ip,el) - np = material_phaseAt(1,neighbor_el) - no = material_phasememberAt(1,neighbor_ip,neighbor_el) - - opposite_neighbor = n + mod(n,2) - mod(n+1,2) - opposite_el = IPneighborhood(1,opposite_neighbor,ip,el) - opposite_ip = IPneighborhood(2,opposite_neighbor,ip,el) - opposite_n = IPneighborhood(3,opposite_neighbor,ip,el) - - if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient - neighbor_instance = phase_plasticityInstance(material_phaseAt(1,neighbor_el)) - neighbor_F = F(1:3,1:3,1,neighbor_ip,neighbor_el) - neighbor_Fe = matmul(neighbor_F, math_inv33(Fp(1:3,1:3,1,neighbor_ip,neighbor_el))) - Favg = 0.5_pReal * (my_F + neighbor_F) - else ! if no neighbor, take my value as average - Favg = my_F - endif - - neighbor_v0 = 0.0_pReal ! needed for check of sign change in flux density below - - !* FLUX FROM MY NEIGHBOR TO ME - !* This is only considered, if I have a neighbor of nonlocal plasticity - !* (also nonlocal constitutive law with local properties) that is at least a little bit - !* compatible. - !* If it's not at all compatible, no flux is arriving, because everything is dammed in front of - !* my neighbor's interface. - !* The entering flux from my neighbor will be distributed on my slip systems according to the - !* compatibility - if (neighbor_n > 0) then - if (phase_plasticity(material_phaseAt(1,neighbor_el)) == PLASTICITY_NONLOCAL_ID .and. & - any(compatibility(:,:,:,n,ip,el) > 0.0_pReal)) then - - forall (s = 1:ns, t = 1:4) - neighbor_v0(s,t) = plasticState(np)%state0(iV (s,t,neighbor_instance),no) - neighbor_rhoSgl0(s,t) = max(plasticState(np)%state0(iRhoU(s,t,neighbor_instance),no),0.0_pReal) - endforall - - where (neighbor_rhoSgl0 * IPvolume(neighbor_ip,neighbor_el) ** 0.667_pReal < prm%significantN & - .or. neighbor_rhoSgl0 < prm%significantRho) & - neighbor_rhoSgl0 = 0.0_pReal - normal_neighbor2me_defConf = math_det33(Favg) * matmul(math_inv33(transpose(Favg)), & - IPareaNormal(1:3,neighbor_n,neighbor_ip,neighbor_el)) ! normal of the interface in (average) deformed configuration (pointing neighbor => me) - normal_neighbor2me = matmul(transpose(neighbor_Fe), normal_neighbor2me_defConf) & - / math_det33(neighbor_Fe) ! interface normal in the lattice configuration of my neighbor - area = IParea(neighbor_n,neighbor_ip,neighbor_el) * norm2(normal_neighbor2me) - normal_neighbor2me = normal_neighbor2me / norm2(normal_neighbor2me) ! normalize the surface normal to unit length - do s = 1,ns - do t = 1,4 - c = (t + 1) / 2 - topp = t + mod(t,2) - mod(t+1,2) - if (neighbor_v0(s,t) * math_inner(m(1:3,s,t), normal_neighbor2me) > 0.0_pReal & ! flux from my neighbor to me == entering flux for me - .and. v0(s,t) * neighbor_v0(s,t) >= 0.0_pReal ) then ! ... only if no sign change in flux density - lineLength = neighbor_rhoSgl0(s,t) * neighbor_v0(s,t) & - * math_inner(m(1:3,s,t), normal_neighbor2me) * area ! positive line length that wants to enter through this interface - where (compatibility(c,:,s,n,ip,el) > 0.0_pReal) & - rhoDotFlux(:,t) = rhoDotFlux(1:ns,t) & - + lineLength/IPvolume(ip,el)*compatibility(c,:,s,n,ip,el)**2.0_pReal ! transferring to equally signed mobile dislocation type - where (compatibility(c,:,s,n,ip,el) < 0.0_pReal) & - rhoDotFlux(:,topp) = rhoDotFlux(:,topp) & - + lineLength/IPvolume(ip,el)*compatibility(c,:,s,n,ip,el)**2.0_pReal ! transferring to opposite signed mobile dislocation type - - endif - enddo - enddo - endif; endif - - - !* FLUX FROM ME TO MY NEIGHBOR - !* This is not considered, if my opposite neighbor has a different constitutive law than nonlocal (still considered for nonlocal law with local properties). - !* Then, we assume, that the opposite(!) neighbor sends an equal amount of dislocations to me. - !* So the net flux in the direction of my neighbor is equal to zero: - !* leaving flux to neighbor == entering flux from opposite neighbor - !* In case of reduced transmissivity, part of the leaving flux is stored as dead dislocation density. - !* That means for an interface of zero transmissivity the leaving flux is fully converted to dead dislocations. - if (opposite_n > 0) then - if (phase_plasticity(material_phaseAt(1,opposite_el)) == PLASTICITY_NONLOCAL_ID) then - - normal_me2neighbor_defConf = math_det33(Favg) & - * matmul(math_inv33(transpose(Favg)),IPareaNormal(1:3,n,ip,el)) ! normal of the interface in (average) deformed configuration (pointing me => neighbor) - normal_me2neighbor = matmul(transpose(my_Fe), normal_me2neighbor_defConf) & - / math_det33(my_Fe) ! interface normal in my lattice configuration - area = IParea(n,ip,el) * norm2(normal_me2neighbor) - normal_me2neighbor = normal_me2neighbor / norm2(normal_me2neighbor) ! normalize the surface normal to unit length - do s = 1,ns - do t = 1,4 - c = (t + 1) / 2 - if (v0(s,t) * math_inner(m(1:3,s,t), normal_me2neighbor) > 0.0_pReal ) then ! flux from me to my neighbor == leaving flux for me (might also be a pure flux from my mobile density to dead density if interface not at all transmissive) - if (v0(s,t) * neighbor_v0(s,t) >= 0.0_pReal) then ! no sign change in flux density - transmissivity = sum(compatibility(c,:,s,n,ip,el)**2.0_pReal) ! overall transmissivity from this slip system to my neighbor - else ! sign change in flux density means sign change in stress which does not allow for dislocations to arive at the neighbor - transmissivity = 0.0_pReal - endif - lineLength = my_rhoSgl0(s,t) * v0(s,t) & - * math_inner(m(1:3,s,t), normal_me2neighbor) * area ! positive line length of mobiles that wants to leave through this interface - rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / IPvolume(ip,el) ! subtract dislocation flux from current type - rhoDotFlux(s,t+4) = rhoDotFlux(s,t+4) & - + lineLength / IPvolume(ip,el) * (1.0_pReal - transmissivity) & - * sign(1.0_pReal, v0(s,t)) ! dislocation flux that is not able to leave through interface (because of low transmissivity) will remain as immobile single density at the material point - endif - enddo - enddo - endif; endif - - enddo neighbors - endif - - !**************************************************************************** !*** calculate dipole formation and annihilation @@ -1454,7 +1276,7 @@ function plastic_nonlocal_dotState2(F,Fp,timestep, instance,of,ip,el) result(rh write(6,'(a)') '<< CONST >> enforcing cutback !!!' endif #endif - plasticState(ph)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN) ! -> return NaN and, hence, enforce cutback + rhoDotFlux = IEEE_value(1.0_pReal,IEEE_quiet_NaN) ! enforce cutback return endif From 520a484df2d81fd3c24e443809c4453082afe67c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 3 May 2020 13:09:49 +0200 Subject: [PATCH 5/9] reasonable name --- src/constitutive_plastic_nonlocal.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index 08d14a5ce..f1994444e 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -1117,7 +1117,7 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, & - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) & - rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have - rhoDot = plastic_nonlocal_dotState2(F,Fp,timestep, instance,of,ip,el) & + rhoDot = rhoDotFlux(F,Fp,timestep, instance,of,ip,el) & + rhoDotMultiplication & + rhoDotSingle2DipoleGlide & + rhoDotAthermalAnnihilation & @@ -1169,7 +1169,7 @@ end subroutine plastic_nonlocal_dotState !--------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !--------------------------------------------------------------------------------------------------- -function plastic_nonlocal_dotState2(F,Fp,timestep, instance,of,ip,el) result(rhoDotFlux) +function rhoDotFlux(F,Fp,timestep, instance,of,ip,el) real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: & F, & !< elastic deformation gradient @@ -1404,7 +1404,7 @@ function plastic_nonlocal_dotState2(F,Fp,timestep, instance,of,ip,el) result(rh end associate -end function plastic_nonlocal_dotState2 +end function rhoDotFlux !-------------------------------------------------------------------------------------------------- From 76d37b5bcd7eba064f9f1f19d36eb9fea37fd9c3 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 3 May 2020 20:38:25 +0200 Subject: [PATCH 6/9] flux debug would need temporary variable --- src/constitutive_plastic_nonlocal.f90 | 23 ----------------------- 1 file changed, 23 deletions(-) diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index f1994444e..cc5655e78 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -1123,29 +1123,6 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, & + rhoDotAthermalAnnihilation & + rhoDotThermalAnnihilation -#ifdef DEBUG - if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0 & - .and. ((debug_e == el .and. debug_i == ip)& - .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0 )) then - write(6,'(a,/,4(12x,12(e12.5,1x),/))') '<< CONST >> dislocation multiplication', & - rhoDotMultiplication(:,1:4) * timestep - write(6,'(a,/,8(12x,12(e12.5,1x),/))') '<< CONST >> dislocation flux', & - rhoDotFlux(:,1:8) * timestep - write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> dipole formation by glide', & - rhoDotSingle2DipoleGlide * timestep - write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> athermal dipole annihilation', & - rhoDotAthermalAnnihilation * timestep - write(6,'(a,/,2(12x,12(e12.5,1x),/))') '<< CONST >> thermally activated dipole annihilation', & - rhoDotThermalAnnihilation(:,9:10) * timestep - write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> total density change', & - rhoDot * timestep - write(6,'(a,/,10(12x,12(f12.5,1x),/))') '<< CONST >> relative density change', & - rhoDot(:,1:8) * timestep / (abs(stt%rho(:,sgl))+1.0e-10), & - rhoDot(:,9:10) * timestep / (stt%rho(:,dip)+1.0e-10) - write(6,*) - endif -#endif - if ( any(rho(:,mob) + rhoDot(:,1:4) * timestep < -prm%atol_rho) & .or. any(rho(:,dip) + rhoDot(:,9:10) * timestep < -prm%atol_rho)) then From 4dcec8b309d11a61d2c8fef77c75050b2ebef0ad Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 5 May 2020 11:18:58 +0200 Subject: [PATCH 7/9] not needed --- src/constitutive_plastic_nonlocal.f90 | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index cc5655e78..fefb668eb 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -1209,10 +1209,6 @@ function rhoDotFlux(F,Fp,timestep, instance,of,ip,el) lineLength !< dislocation line length leaving the current interface ph = material_phaseAt(1,el) - if (timestep <= 0.0_pReal) then - plasticState(ph)%dotState = 0.0_pReal - return - endif associate(prm => param(instance), & dst => microstructure(instance), & From a5b78dc30af2211b1188407b3b1e277f1c0a0855 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 5 May 2020 11:20:16 +0200 Subject: [PATCH 8/9] potential glitch --- src/constitutive_plastic_nonlocal.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index fefb668eb..65fea24b2 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -1223,7 +1223,7 @@ function rhoDotFlux(F,Fp,timestep, instance,of,ip,el) rho0 = getRho0(instance,of,ip,el) my_rhoSgl0 = rho0(:,sgl) - forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,instance),of) + forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,instance),of) !ToDo: MD: I think we should use state0 here gdot = rhoSgl(:,1:4) * v * spread(prm%burgers,2,4) From 58537c478d3edc95f26d392b7ae3cfb77ecdce00 Mon Sep 17 00:00:00 2001 From: Test User Date: Wed, 6 May 2020 16:01:35 +0200 Subject: [PATCH 9/9] [skip ci] updated version information after successful test of v2.0.3-2412-g0d03c469 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 6f8955f5c..f955915ff 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-2402-gb88f5ec0 +v2.0.3-2412-g0d03c469