From 97db70cf232f8a596c843d7ba87f6ba361b10f43 Mon Sep 17 00:00:00 2001 From: Christoph Kords Date: Wed, 17 Feb 2010 13:21:36 +0000 Subject: [PATCH] constitutive_nonlocal: - reworked contribution of immobile dislocation density for rate equations - flux is now calculated on the basis of interpolated velocities and densities at the interface; both incoming and outgoing fluxes are considered, so every material point only changes his own dotState - dislocation velocity is now globally defined and calculated by subroutine constitutive_nonlocal_kinetics; the subroutine is called inside _LpAndItsTangent as well as _microstructure; therefore, microstructure now needs Tstar_v as additional input; in the future one should perhaps create a subroutine constitutive_kinetics that calls constitutive_nonlocal_kinetics separately, to clearly distinguish between microstructural and kinetic variables - better use flux density vector as output variable instead of scalar flux values for each interface - added output variables internal and external resolved stress crystallite: - added flag to force local stiffness calculation in case of nonlocal model - misorientation angle is explicitly set to zero when no neighbor can be found debug: - added flag "selectiveDebugger" that is used when debugging statements should only affect a specific element, ip and grain; these are specified with the new variables debug_e, debug_i and debug_g - debugger can now be used in its original sense --- code/CPFEM.f90 | 13 +- code/constitutive.f90 | 9 +- code/constitutive_nonlocal.f90 | 695 +++++++++++++++++---------------- code/crystallite.f90 | 74 ++-- code/debug.f90 | 4 + code/homogenization.f90 | 9 +- code/material.config | 43 +- 7 files changed, 432 insertions(+), 415 deletions(-) diff --git a/code/CPFEM.f90 b/code/CPFEM.f90 index 5d9e5fb42..d31b4deca 100644 --- a/code/CPFEM.f90 +++ b/code/CPFEM.f90 @@ -73,7 +73,10 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt use numerics, only: numerics_init, & relevantStrain, & iJacoStiffness - use debug, only: debug_init + use debug, only: debug_init, & + debug_g, & + debug_i, & + debug_e use FEsolving, only: FE_init, & parallelExecution, & outdatedFFN1, & @@ -81,6 +84,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt cycleCounter, & theInc, & theTime, & + theDelta, & FEsolving_execElem, & FEsolving_execIP use math, only: math_init, & @@ -205,8 +209,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt !$OMP CRITICAL (write2out) write(6,*) write(6,*) '#####################################' - write(6,'(a10,x,f8.4,x,a10,x,i4,x,a10,x,i3,x,a16,x,i2,x,a16,x,i2)') & - 'theTime',theTime,'theInc',theInc,'cycleCounter',cycleCounter,'computationMode',mode + write(6,'(a10,x,f8.4,x,a10,x,f8.4,x,a10,x,i6,x,a10,x,i3,x,a16,x,i2,x,a16,x,i2)') & + 'theTime',theTime,'theDelta',theDelta,'theInc',theInc,'cycleCounter',cycleCounter,'computationMode',mode write(6,*) '#####################################' call flush (6) !$OMP END CRITICAL (write2out) @@ -228,7 +232,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt k = 1:mesh_NcpElems ) & constitutive_state0(i,j,k)%p = constitutive_state(i,j,k)%p ! microstructure of crystallites !$OMP CRITICAL (write2out) - write(6,'(a10,/,4(3(e20.8,x),/))') 'aged state',constitutive_state(1,1,1)%p + write(6,'(a,3(x,i4),/,4(3(e20.8,x),/))') 'aged state at', debug_g, debug_i, debug_e, & + constitutive_state(debug_g,debug_i,debug_e)%p !$OMP END CRITICAL (write2out) do k = 1,mesh_NcpElems do j = 1,mesh_maxNips diff --git a/code/constitutive.f90 b/code/constitutive.f90 index 2caf68665..412469a94 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -45,7 +45,7 @@ subroutine constitutive_init() !* Module initialization * !************************************** use prec, only: pReal,pInt - use debug, only: debugger + use debug, only: debugger, selectiveDebugger, debug_e, debug_i, debug_g use IO, only: IO_error, IO_open_file, IO_open_jobFile use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips use material @@ -125,7 +125,7 @@ subroutine constitutive_init() myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = 1,FE_Nips(mesh_element(2,e)) ! loop over IPs do g = 1,myNgrains ! loop over grains - debugger = (e == 1 .and. i == 1 .and. g == 1) + selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) myInstance = phase_constitutionInstance(material_phase(g,i,e)) select case(phase_constitution(material_phase(g,i,e))) @@ -269,7 +269,7 @@ return endfunction -subroutine constitutive_microstructure(Temperature,Fe,Fp,ipc,ip,el) +subroutine constitutive_microstructure(Temperature,Tstar_v,Fe,Fp,ipc,ip,el) !********************************************************************* !* This function calculates from state needed variables * !* INPUT: * @@ -294,6 +294,7 @@ subroutine constitutive_microstructure(Temperature,Fe,Fp,ipc,ip,el) !* Definition of variables integer(pInt), intent(in) :: ipc,ip,el real(pReal), intent(in) :: Temperature +real(pReal), dimension(6) :: Tstar_v real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: Fe, Fp select case (phase_constitution(material_phase(ipc,ip,el))) @@ -308,7 +309,7 @@ real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) call constitutive_dislotwin_microstructure(Temperature,constitutive_state,ipc,ip,el) case (constitutive_nonlocal_label) - call constitutive_nonlocal_microstructure(constitutive_state, Temperature, Fe, Fp, ipc, ip, el) + call constitutive_nonlocal_microstructure(constitutive_state, Temperature, Tstar_v, Fe, Fp, ipc, ip, el) end select diff --git a/code/constitutive_nonlocal.f90 b/code/constitutive_nonlocal.f90 index c1b5b27e0..d054aa0c6 100644 --- a/code/constitutive_nonlocal.f90 +++ b/code/constitutive_nonlocal.f90 @@ -17,19 +17,19 @@ implicit none !* Definition of parameters character (len=*), parameter :: constitutive_nonlocal_label = 'nonlocal' -character(len=20), dimension(10), parameter :: constitutive_nonlocal_listBasicStates = (/'rhoSglEdgePos ', & - 'rhoSglEdgeNeg ', & - 'rhoSglScrewPos ', & - 'rhoSglScrewNeg ', & - 'rhoSglEdgePosUsed ', & - 'rhoSglEdgeNegUsed ', & - 'rhoSglScrewPosUsed ', & - 'rhoSglScrewNegUsed ', & - 'rhoDipEdge ', & - 'rhoDipScrew ' /) ! list of "basic" microstructural state variables that are independent from other state variables -character(len=20), dimension(3), parameter :: constitutive_nonlocal_listDependentStates = (/'rhoForest ', & - 'tauSlipThreshold ', & - 'Tdislocation_v ' /) ! list of microstructural state variables that depend on other state variables +character(len=22), dimension(10), parameter :: constitutive_nonlocal_listBasicStates = (/'rhoSglEdgePosMobile ', & + 'rhoSglEdgeNegMobile ', & + 'rhoSglScrewPosMobile ', & + 'rhoSglScrewNegMobile ', & + 'rhoSglEdgePosImmobile ', & + 'rhoSglEdgeNegImmobile ', & + 'rhoSglScrewPosImmobile', & + 'rhoSglScrewNegImmobile', & + 'rhoDipEdge ', & + 'rhoDipScrew ' /) ! list of "basic" microstructural state variables that are independent from other state variables +character(len=15), dimension(3), parameter :: constitutive_nonlocal_listDependentStates = (/'rhoForest ', & + 'tauThreshold ', & + 'Tdislocation_v ' /) ! list of microstructural state variables that depend on other state variables real(pReal), parameter :: kB = 1.38e-23_pReal ! Physical parameter, Boltzmann constant in J/Kelvin !* Definition of global variables @@ -78,7 +78,7 @@ real(pReal), dimension(:,:), allocatable :: constitutive_nonlocal_ constitutive_nonlocal_dLowerScrewPerSlipFamily, & ! minimum stable screw dipole height for each family and instance constitutive_nonlocal_dLowerScrewPerSlipSystem, & ! minimum stable screw dipole height for each slip system and instance constitutive_nonlocal_interactionSlipSlip ! coefficients for slip-slip interaction for each interaction type and instance - +real(pReal), dimension(:,:,:,:,:), allocatable :: constitutive_nonlocal_v ! dislocation velocity real(pReal), dimension(:,:,:), allocatable :: constitutive_nonlocal_forestProjectionEdge, & ! matrix of forest projections of edge dislocations for each instance constitutive_nonlocal_forestProjectionScrew, & ! matrix of forest projections of screw dislocations for each instance constitutive_nonlocal_interactionMatrixSlipSlip ! interaction matrix of the different slip systems for each instance @@ -114,7 +114,10 @@ use IO, only: IO_lc, & IO_floatValue, & IO_intValue, & IO_error -use material, only: phase_constitution, & +use mesh, only: mesh_NcpElems, & + mesh_maxNips +use material, only: homogenization_maxNgrains, & + phase_constitution, & phase_constitutionInstance, & phase_Noutput use lattice, only: lattice_maxNslipFamily, & @@ -408,6 +411,8 @@ constitutive_nonlocal_forestProjectionScrew = 0.0_pReal allocate(constitutive_nonlocal_interactionMatrixSlipSlip(maxTotalNslip, maxTotalNslip, maxNinstance)) constitutive_nonlocal_interactionMatrixSlipSlip = 0.0_pReal +allocate(constitutive_nonlocal_v(maxTotalNslip, 4, homogenization_maxNgrains, mesh_maxNips, mesh_NcpElems)) +constitutive_nonlocal_v = 0.0_pReal do i = 1,maxNinstance @@ -472,6 +477,8 @@ do i = 1,maxNinstance 'rho_forest', & 'shearrate', & 'resolvedstress', & + 'resolvedstress_internal', & + 'resolvedstress_external', & 'resistance', & 'rho_dot', & 'rho_dot_sgl', & @@ -483,35 +490,18 @@ do i = 1,maxNinstance 'rho_dot_dip2sgl', & 'rho_dot_ann_ath', & 'rho_dot_ann_the', & - 'rho_dot_flux', & - 'rho_dot_flux_ep', & - 'rho_dot_flux_en', & - 'rho_dot_flux_sp', & - 'rho_dot_flux_sn', & - 'rho_dot_flux_n1_ep', & - 'rho_dot_flux_n1_en', & - 'rho_dot_flux_n1_sp', & - 'rho_dot_flux_n1_sn', & - 'rho_dot_flux_n2_ep', & - 'rho_dot_flux_n2_en', & - 'rho_dot_flux_n2_sp', & - 'rho_dot_flux_n2_sn', & - 'rho_dot_flux_n3_ep', & - 'rho_dot_flux_n3_en', & - 'rho_dot_flux_n3_sp', & - 'rho_dot_flux_n3_sn', & - 'rho_dot_flux_n4_ep', & - 'rho_dot_flux_n4_en', & - 'rho_dot_flux_n4_sp', & - 'rho_dot_flux_n4_sn', & - 'rho_dot_flux_n5_ep', & - 'rho_dot_flux_n5_en', & - 'rho_dot_flux_n5_sp', & - 'rho_dot_flux_n5_sn', & - 'rho_dot_flux_n6_ep', & - 'rho_dot_flux_n6_en', & - 'rho_dot_flux_n6_sp', & - 'rho_dot_flux_n6_sn', & + 'fluxdensity_edge_pos_x', & + 'fluxdensity_edge_pos_y', & + 'fluxdensity_edge_pos_z', & + 'fluxdensity_edge_neg_x', & + 'fluxdensity_edge_neg_y', & + 'fluxdensity_edge_neg_z', & + 'fluxdensity_screw_pos_x', & + 'fluxdensity_screw_pos_y', & + 'fluxdensity_screw_pos_z', & + 'fluxdensity_screw_neg_x', & + 'fluxdensity_screw_neg_y', & + 'fluxdensity_screw_neg_z', & 'd_upper_edge', & 'd_upper_screw', & 'd_upper_dot_edge', & @@ -737,7 +727,7 @@ endfunction !********************************************************************* !* calculates quantities characterizing the microstructure * !********************************************************************* -subroutine constitutive_nonlocal_microstructure(state, Temperature, Fe, Fp, g, ip, el) +subroutine constitutive_nonlocal_microstructure(state, Temperature, Tstar_v, Fe, Fp, g, ip, el) use prec, only: pReal, & pInt, & @@ -751,7 +741,8 @@ use math, only: math_Plain3333to99, & math_inv3x3, & math_det3x3, & pi -use debug, only: debugger +use debug, only: debugger, & + selectiveDebugger use mesh, only: mesh_NcpElems, & mesh_maxNips, & mesh_element, & @@ -781,6 +772,9 @@ real(pReal), intent(in) :: Temperature ! temperature real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & Fe, & ! elastic deformation gradient Fp ! plastic deformation gradient +real(pReal), dimension(6), intent(in) :: & + Tstar_v ! 2nd Piola-Kirchhoff stress in Mandel notation + !*** input/output variables type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(inout) :: & @@ -822,7 +816,8 @@ real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstan rhoDip ! dipole dislocation density (edge, screw) real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: & rhoForest, & ! forest dislocation density - tauSlipThreshold, & ! threshold shear stress + tauThreshold, & ! threshold shear stress + tau, & ! resolved shear stress neighboring_rhoEdgeExcess, &! edge excess dislocation density of my neighbor neighboring_rhoScrewExcess,&! screw excess dislocation density of my neighbor neighboring_Nedge, & ! total number of edge excess dislocations in my neighbor @@ -857,10 +852,10 @@ forall (s = 1:ns) & !*** calculate the threshold shear stress for dislocation slip forall (s = 1:ns) & - tauSlipThreshold(s) = constitutive_nonlocal_Gmod(myInstance) & - * constitutive_nonlocal_burgersPerSlipSystem(s, myInstance) & - * sqrt( dot_product( ( sum(abs(rhoSgl),2) + sum(abs(rhoDip),2) ), & - constitutive_nonlocal_interactionMatrixSlipSlip(s, 1:ns, myInstance) ) ) + tauThreshold(s) = constitutive_nonlocal_Gmod(myInstance) & + * constitutive_nonlocal_burgersPerSlipSystem(s, myInstance) & + * sqrt( dot_product( ( sum(abs(rhoSgl),2) + sum(abs(rhoDip),2) ), & + constitutive_nonlocal_interactionMatrixSlipSlip(s, 1:ns, myInstance) ) ) ! if (debugger) write(6,'(a26,3(i3,x),/,12(f10.5,x),/)') 'tauSlipThreshold / MPa at ',g,ip,el, tauSlipThreshold/1e6 @@ -932,14 +927,103 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el)) enddo enddo - !********************************************************************** !*** set dependent states state(g,ip,el)%p(10*ns+1:11*ns) = rhoForest -state(g,ip,el)%p(11*ns+1:12*ns) = tauSlipThreshold +state(g,ip,el)%p(11*ns+1:12*ns) = tauThreshold state(g,ip,el)%p(12*ns+1:12*ns+6) = Tdislocation_v + +!********************************************************************** +!*** calculate the dislocation velocity + +call constitutive_nonlocal_kinetics(Tstar_v, Temperature, state, g, ip, el) + +endsubroutine + + + +!********************************************************************* +!* calculates kinetics * +!********************************************************************* +subroutine constitutive_nonlocal_kinetics(Tstar_v, Temperature, state, g, ip, el) + +use prec, only: pReal, & + pInt, & + p_vec +use math, only: math_mul6x6, & + math_Mandel6to33 +use debug, only: debugger, & + selectiveDebugger +use mesh, only: mesh_NcpElems, & + mesh_maxNips +use material, only: homogenization_maxNgrains, & + material_phase, & + phase_constitutionInstance +use lattice, only: lattice_Sslip, & + lattice_Sslip_v + +implicit none + +!*** input variables +integer(pInt), intent(in) :: g, & ! current grain number + ip, & ! current integration point + el ! current element number +real(pReal), intent(in) :: Temperature ! temperature +type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + state ! microstructural state +real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola-Kirchhoff stress in Mandel notation + +!*** output variables + +!*** local variables +integer(pInt) myInstance, & ! current instance of this constitution + myStructure, & ! current lattice structure + ns, & ! short notation for the total number of active slip systems + t, & ! dislocation type + s ! index of my current slip system +real(pReal), dimension(6) :: Tdislocation_v ! dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress +real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: & + tauThreshold, & ! threshold shear stress + tau ! resolved shear stress + + +myInstance = phase_constitutionInstance(material_phase(g,ip,el)) +myStructure = constitutive_nonlocal_structure(myInstance) +ns = constitutive_nonlocal_totalNslip(myInstance) + +tauThreshold = state(g,ip,el)%p(11*ns+1:12*ns) +Tdislocation_v = state(g,ip,el)%p(12*ns+1:12*ns+6) + +tau = 0.0_pReal +constitutive_nonlocal_v(:,:,g,ip,el) = 0.0_pReal + +do s = 1,ns + if ((tauThreshold(s) > 0.0_pReal) .and. (Temperature > 0.0_pReal)) then + + tau(s) = math_mul6x6( Tstar_v + Tdislocation_v, & + lattice_Sslip_v(:,constitutive_nonlocal_slipSystemLattice(s,myInstance),myStructure) ) + + constitutive_nonlocal_v(s,:,g,ip,el) = constitutive_nonlocal_v0PerSlipSystem(s,myInstance) & + * exp( - constitutive_nonlocal_Q0(myInstance) / ( kB * Temperature) * (1.0_pReal - (abs(tau(s))/tauThreshold(s)) ) ) & + * sign(1.0_pReal,tau(s)) + + endif +enddo + +if (selectiveDebugger) then + !$OMP CRITICAL (write2out) + write(6,*) '::: kinetics',g,ip,el + write(6,*) + write(6,'(a,/,3(3(f12.3,x)/))') 'Tdislocation / MPa', math_Mandel6to33(Tdislocation_v/1e6) + write(6,'(a,/,3(3(f12.3,x)/))') 'Tstar / MPa', math_Mandel6to33(Tstar_v/1e6) + write(6,'(a,/,12(f12.5,x),/)') 'tau / MPa', tau/1e6_pReal + write(6,'(a,/,12(f12.5,x),/)') 'tauThreshold / MPa', tauThreshold/1e6_pReal + write(6,'(a,/,4(12(f12.5,x),/))') 'v', constitutive_nonlocal_v(:,:,g,ip,el) + !$OMPEND CRITICAL (write2out) +endif + endsubroutine @@ -955,7 +1039,8 @@ use prec, only: pReal, & use math, only: math_Plain3333to99, & math_mul6x6, & math_Mandel6to33 -use debug, only: debugger +use debug, only: debugger, & + selectiveDebugger use mesh, only: mesh_NcpElems, & mesh_maxNips use material, only: homogenization_maxNgrains, & @@ -990,24 +1075,20 @@ integer(pInt) myInstance, & ! curren t, & ! dislocation type s, & ! index of my current slip system sLattice ! index of my current slip system according to lattice order -real(pReal), dimension(6) :: Tdislocation_v ! dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress real(pReal), dimension(3,3,3,3) :: dLp_dTstar3333 ! derivative of Lp with respect to Tstar (3x3x3x3 matrix) real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),8) :: & rhoSgl ! single dislocation densities (including used) +real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4) :: & + gdot ! shear rate per dislocation type real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: & - rhoForest, & ! forest dislocation density - tauSlipThreshold, & ! threshold shear stress - tauSlip, & ! resolved shear stress - gdotSlip, & ! shear rate - dgdot_dtauSlip, & ! derivative of the shear rate with respect to the shear stress - v ! dislocation velocity + tauThreshold, & ! threshold shear stress + gdotTotal, & ! shear rate + dgdotTotal_dtau ! derivative of the shear rate with respect to the shear stress !*** initialize local variables -v = 0.0_pReal -tauSlip = 0.0_pReal -gdotSlip = 0.0_pReal +gdot = 0.0_pReal Lp = 0.0_pReal dLp_dTstar3333 = 0.0_pReal @@ -1017,54 +1098,45 @@ ns = constitutive_nonlocal_totalNslip(myInstance) !*** shortcut to state variables -forall (t = 1:8) rhoSgl(:,t) = state(g,ip,el)%p((t-1)*ns+1:t*ns) -rhoForest = state(g,ip,el)%p(10*ns+1:11*ns) -tauSlipThreshold = state(g,ip,el)%p(11*ns+1:12*ns) -Tdislocation_v = state(g,ip,el)%p(12*ns+1:12*ns+6) +forall (t = 1:8) & + rhoSgl(:,t) = state(g,ip,el)%p((t-1)*ns+1:t*ns) +forall (s = 1:ns, t = 5:8, rhoSgl(s,t) * constitutive_nonlocal_v(s,t-4,g,ip,el) < 0.0_pReal) & ! contribution of used rho for changing sign of v + rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t)) -!*** calculation of resolved stress +tauThreshold = state(g,ip,el)%p(11*ns+1:12*ns) -forall (s = 1:ns) & - tauSlip(s) = math_mul6x6( Tstar_v + Tdislocation_v, & - lattice_Sslip_v(:,constitutive_nonlocal_slipSystemLattice(s,myInstance),myStructure) ) +call constitutive_nonlocal_kinetics(Tstar_v, Temperature, state, g, ip, el) ! update dislocation velocity !*** Calculation of gdot and its tangent -v = constitutive_nonlocal_v0PerSlipSystem(:,myInstance) & - * exp( - constitutive_nonlocal_Q0(myInstance) / ( kB * Temperature) * (1.0_pReal - (abs(tauSlip)/tauSlipThreshold) ) ) & - * sign(1.0_pReal,tauSlip) +forall (t = 1:4 ) & + gdot(:,t) = rhoSgl(:,t) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) * constitutive_nonlocal_v(:,t,g,ip,el) +gdotTotal = sum(gdot,2) -gdotSlip = sum(rhoSgl(:,1:4),2) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) * v -forall (s = 1:ns, t = 5:8, rhoSgl(s,t) * v(s) < 0.0_pReal) & ! contribution of used rho for changing sign of v - gdotSlip(s) = gdotSlip(s) + abs(rhoSgl(s,t)) * constitutive_nonlocal_burgersPerSlipSystem(s,myInstance) * v(s) - -dgdot_dtauSlip = abs(gdotSlip) * constitutive_nonlocal_Q0(myInstance) / ( kB * Temperature * tauSlipThreshold ) +dgdotTotal_dtau = abs(gdotTotal) * constitutive_nonlocal_Q0(myInstance) / ( kB * Temperature * tauThreshold ) !*** Calculation of Lp and its tangent do s = 1,ns sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance) - Lp = Lp + gdotSlip(s) * lattice_Sslip(:,:,sLattice,myStructure) + Lp = Lp + gdotTotal(s) * lattice_Sslip(:,:,sLattice,myStructure) forall (i=1:3,j=1:3,k=1:3,l=1:3) & - dLp_dTstar3333(i,j,k,l) = dLp_dTstar3333(i,j,k,l) + dgdot_dtauSlip(s) * lattice_Sslip(i,j, sLattice,myStructure) & - * lattice_Sslip(k,l, sLattice,myStructure) + dLp_dTstar3333(i,j,k,l) = dLp_dTstar3333(i,j,k,l) + dgdotTotal_dtau(s) * lattice_Sslip(i,j, sLattice,myStructure) & + * lattice_Sslip(k,l, sLattice,myStructure) enddo dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333) -if (debugger) then +if (selectiveDebugger) then !$OMP CRITICAL (write2out) write(6,*) '::: LpandItsTangent',g,ip,el write(6,*) - write(6,'(a,/,3(3(f12.3,x)/))') 'Tdislocation / MPa', math_Mandel6to33(Tdislocation_v/1e6) - write(6,'(a,/,3(3(f12.3,x)/))') 'Tstar / MPa', math_Mandel6to33(Tstar_v/1e6) - write(6,'(a,/,12(f12.5,x),/)') 'tauSlip / MPa', tauSlip/1e6_pReal - ! write(6,'(a,/,12(f12.5,x),/)') 'tauSlipThreshold / MPa', tauSlipThreshold/1e6_pReal - ! write(6,'(a,/,12(f12.5,x),/)') 'v', v - ! write(6,'(a,/,12(f12.5,x),/)') 'gdot total /1e-3',gdotSlip*1e3_pReal - ! write(6,'(a,/,3(3(f12.7,x)/))') 'Lp',Lp + ! write(6,'(a,/,12(f12.5,x),/)') 'v', constitutive_nonlocal_v(:,t,g,ip,el) + write(6,'(a,/,12(f12.5,x),/)') 'gdot /1e-3',gdot*1e3_pReal + write(6,'(a,/,12(f12.5,x),/)') 'gdot total /1e-3',gdotTotal*1e3_pReal + write(6,'(a,/,3(3(f12.7,x)/))') 'Lp',Lp ! call flush(6) !$OMPEND CRITICAL (write2out) endif @@ -1082,7 +1154,8 @@ subroutine constitutive_nonlocal_dotState(dotState, Tstar_v, previousTstar_v, Fe use prec, only: pReal, & pInt, & p_vec -use debug, only: debugger +use debug, only: debugger, & + selectiveDebugger use math, only: math_norm3, & math_mul6x6, & math_mul3x3, & @@ -1154,14 +1227,12 @@ real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstan totalRhoDotSgl, & ! total rate of change of single dislocation densities thisRhoDotSgl ! rate of change of single dislocation densities for this mechanism real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4) :: & - gdot, & ! shear rates - lineLength ! dislocation line length leaving the current interface + gdot ! shear rates real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: & rhoForest, & ! forest dislocation density - tauSlipThreshold, & ! threshold shear stress - tauSlip, & ! current resolved shear stress - previousTauSlip, & ! previous resolved shear stress - v, & ! dislocation velocity + tauThreshold, & ! threshold shear stress + tau, & ! current resolved shear stress + previousTau, & ! previous resolved shear stress invLambda, & ! inverse of mean free path for dislocations vClimb ! climb velocity of edge dipoles real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),2) :: & @@ -1184,6 +1255,10 @@ real(pReal), dimension(3) :: surfaceNormal, & ! surface surfaceNormal_currentconf ! surface normal in current configuration real(pReal) area, & ! area of the current interface detFe, & ! determinant of elastic defornmation gradient + transmissivity, & ! transmissivity of interfaces for dislocation flux + rhoAvg, & ! average mobile single dislocation density at interface + vAvg, & ! average dislocation velocity at interface + lineLength, & ! dislocation line length leaving the current interface D ! self diffusion @@ -1191,9 +1266,8 @@ myInstance = phase_constitutionInstance(material_phase(g,ip,el)) myStructure = constitutive_nonlocal_structure(myInstance) ns = constitutive_nonlocal_totalNslip(myInstance) -tauSlip = 0.0_pReal -previousTauSlip = 0.0_pReal -v = 0.0_pReal +tau = 0.0_pReal +previousTau = 0.0_pReal gdot = 0.0_pReal dLower = 0.0_pReal dUpper = 0.0_pReal @@ -1211,7 +1285,7 @@ forall (t = 1:8) previousRhoSgl(:,t) = previousState(g,ip,el)%p((t-1)*ns+1:t*ns) forall (c = 1:2) rhoDip(:,c) = state(g,ip,el)%p((7+c)*ns+1:(8+c)*ns) forall (c = 1:2) previousRhoDip(:,c) = previousState(g,ip,el)%p((7+c)*ns+1:(8+c)*ns) rhoForest = state(g,ip,el)%p(10*ns+1:11*ns) -tauSlipThreshold = state(g,ip,el)%p(11*ns+1:12*ns) +tauThreshold = state(g,ip,el)%p(11*ns+1:12*ns) Tdislocation_v = state(g,ip,el)%p(12*ns+1:12*ns+6) previousTdislocation_v = previousState(g,ip,el)%p(12*ns+1:12*ns+6) @@ -1219,52 +1293,41 @@ previousTdislocation_v = previousState(g,ip,el)%p(12*ns+1:12*ns+6) !**************************************************************************** !*** Calculate shear rate -do s = 1,ns ! loop over slip systems - sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance) - - tauSlip(s) = math_mul6x6( Tstar_v + Tdislocation_v, lattice_Sslip_v(:,sLattice,myStructure) ) - previousTauSlip(s) = math_mul6x6( previousTstar_v + previousTdislocation_v, lattice_Sslip_v(:,sLattice,myStructure) ) - -enddo - -v = constitutive_nonlocal_v0PerSlipSystem(:,myInstance) & - * exp( - constitutive_nonlocal_Q0(myInstance) / ( kB * Temperature) * (1.0_pReal - (abs(tauSlip)/tauSlipThreshold) ) ) & - * sign(1.0_pReal,tauSlip) - forall (t = 1:4) & - gdot(:,t) = rhoSgl(:,t) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) * v -forall (s = 1:ns, t = 1:4, rhoSgl(s,t+4) * v(s) < 0.0_pReal) & ! contribution of used rho for changing sign of v - gdot(s,t) = gdot(s,t) + abs(rhoSgl(s,t+4)) * constitutive_nonlocal_burgersPerSlipSystem(s,myInstance) * v(s) + gdot(:,t) = rhoSgl(:,t) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) * constitutive_nonlocal_v(:,t,g,ip,el) +forall (s = 1:ns, t = 1:4, rhoSgl(s,t+4) * constitutive_nonlocal_v(s,t,g,ip,el) < 0.0_pReal) & ! contribution of used rho for changing sign of v + gdot(s,t) = gdot(s,t) + abs(rhoSgl(s,t+4)) * constitutive_nonlocal_burgersPerSlipSystem(s,myInstance) & + * constitutive_nonlocal_v(s,t,g,ip,el) -if (debugger) then +if (selectiveDebugger) then !$OMP CRITICAL (write2out) write(6,*) '::: constitutive_nonlocal_dotState at ',g,ip,el write(6,*) - write(6,'(a,/,3(3(f12.6,x)/))') 'Tdislocation / MPa', math_Mandel6to33(Tdislocation_v/1e6) - write(6,'(a,/,3(3(f12.6,x)/))') 'Tstar / MPa', math_Mandel6to33(Tstar_v/1e6) - write(6,'(a,/,12(f12.5,x),/)') 'tauSlip / MPa', tauSlip / 1e6_pReal - write(6,'(a,/,12(f12.5,x),/)') 'tauSlipThreshold / MPa', tauSlipThreshold / 1e6_pReal - write(6,'(a,/,12(e12.5,x),/)') 'v', v - write(6,'(a,/,4(12(e12.5,x),/))') 'gdot / 1e-3', gdot*1e3_pReal - ! write(6,'(a,/,(12(f12.5,x),/))') 'gdot total/ 1e-3', sum(gdot,2)*1e3_pReal - ! call flush(6) !$OMPEND CRITICAL (write2out) endif !**************************************************************************** !*** calculate limits for stable dipole height and its rate of change + +do s = 1,ns ! loop over slip systems + sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance) + tau(s) = math_mul6x6( Tstar_v + Tdislocation_v, lattice_Sslip_v(:,sLattice,myStructure) ) + previousTau(s) = math_mul6x6( previousTstar_v + previousTdislocation_v, lattice_Sslip_v(:,sLattice,myStructure) ) + +enddo + dLower(:,1) = constitutive_nonlocal_dLowerEdgePerSlipSystem(:,myInstance) dLower(:,2) = constitutive_nonlocal_dLowerScrewPerSlipSystem(:,myInstance) -dUpper(:,2) = min( constitutive_nonlocal_Gmod(myInstance) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) & - / ( 8.0_pReal * pi * abs(tauSlip) ), & - 1.0_pReal / sqrt( sum(abs(rhoSgl),2)+sum(rhoDip,2) ) ) +dUpper(:,2) = min( 1.0_pReal / sqrt( sum(abs(rhoSgl),2)+sum(rhoDip,2) ), & + constitutive_nonlocal_Gmod(myInstance) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) & + / ( 8.0_pReal * pi * abs(tau) ) ) dUpper(:,1) = dUpper(:,2) / ( 1.0_pReal - constitutive_nonlocal_nu(myInstance) ) -previousDUpper(:,2) = min( constitutive_nonlocal_Gmod(myInstance) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) & - / ( 8.0_pReal * pi * abs(previousTauSlip) ), & - 1.0_pReal / sqrt( sum(abs(previousRhoSgl),2) + sum(previousRhoDip,2) ) ) +previousDUpper(:,2) = min( 1.0_pReal / sqrt( sum(abs(previousRhoSgl),2) + sum(previousRhoDip,2) ), & + constitutive_nonlocal_Gmod(myInstance) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) & + / ( 8.0_pReal * pi * abs(previousTau) ) ) previousDUpper(:,1) = previousDUpper(:,2) / ( 1.0_pReal - constitutive_nonlocal_nu(myInstance) ) if (dt_previous > 0.0_pReal) dUpperDot = (dUpper - previousDUpper) / dt_previous @@ -1273,7 +1336,6 @@ if (debugger) then !$OMP CRITICAL (write2out) ! write(6,'(a,/,2(12(f12.5,x),/))') 'dUpper / micron', dUpper*1e6_pReal ! write(6,'(a,/,2(12(f12.5,x),/))') 'dUpperDot / micron/s', dUpperDot * 1e6_pReal - ! call flush(6) !$OMPEND CRITICAL (write2out) endif @@ -1285,10 +1347,10 @@ thisRhoDotSgl = 0.0_pReal if (timestep > 0.0_pReal) then do t = 1,4 do s = 1,ns - if (rhoSgl(s,t+4) * v(s) < 0.0_pReal) then + if (rhoSgl(s,t+4) * constitutive_nonlocal_v(s,t,g,ip,el) < 0.0_pReal) then thisRhoDotSgl(s,t) = abs(rhoSgl(s,t+4)) / timestep rhoSgl(s,t) = rhoSgl(s,t) + abs(rhoSgl(s,t+4)) - thisRhoDotSgl(s,t+4) = - sign(1.0_pReal,rhoSgl(s,t+4)) * rhoSgl(s,t+4) / timestep + thisRhoDotSgl(s,t+4) = - rhoSgl(s,t+4) / timestep rhoSgl(s,t+4) = 0.0_pReal endif enddo @@ -1297,7 +1359,7 @@ endif totalRhoDotSgl = totalRhoDotSgl + thisRhoDotSgl -if (debugger) then +if (selectiveDebugger) then !$OMP CRITICAL (write2out) write(6,'(a,/,8(12(e12.5,x),/))') 'dislocation remobilization', thisRhoDotSgl * timestep !$OMPEND CRITICAL (write2out) @@ -1319,7 +1381,7 @@ thisRhoDotDip = 0.0_pReal ! dipoles don't multiplicate totalRhoDotSgl = totalRhoDotSgl + thisRhoDotSgl totalRhoDotDip = totalRhoDotDip + thisRhoDotDip -if (debugger) then +if (selectiveDebugger) then !$OMP CRITICAL (write2out) write(6,'(a,/,4(12(e12.5,x),/))') 'dislocation multiplication', thisRhoDotSgl(:,1:4) * timestep !$OMPEND CRITICAL (write2out) @@ -1342,56 +1404,77 @@ detFe = math_det3x3(Fe(:,:,g,ip,el)) do n = 1,FE_NipNeighbors(mesh_element(2,el)) ! loop through my neighbors + transmissivity = constitutive_nonlocal_transmissivity(misorientation(4,n), misorientation(1:3,n)) + neighboring_el = mesh_ipNeighborhood(1,n,ip,el) neighboring_ip = mesh_ipNeighborhood(2,n,ip,el) - ! if neighbor exists, total deformation gradient is averaged over me and my neighbor - if ( neighboring_el > 0 .and. neighboring_ip > 0 ) then + if ( neighboring_el > 0 .and. neighboring_ip > 0 ) then ! if neighbor exists, average deformation gradient neighboring_F = math_mul33x33(Fe(:,:,g,neighboring_ip,neighboring_el), Fp(:,:,g,neighboring_ip,neighboring_el)) Favg = 0.5_pReal * (F + neighboring_F) - else + else ! if no neighbor, take my value as average Favg = F - endif - - ! calculate the normal of the interface in current and lattice configuration - surfaceNormal_currentconf = math_det3x3(Favg) * math_mul33x3(math_inv3x3(transpose(Favg)), mesh_ipAreaNormal(:,n,ip,el)) - surfaceNormal = math_mul33x3(transpose(Fe(:,:,g,ip,el)), surfaceNormal_currentconf) / detFe + endif - ! calculate the area of the interface + surfaceNormal_currentconf = math_det3x3(Favg) * math_mul33x3(math_inv3x3(transpose(Favg)), mesh_ipAreaNormal(:,n,ip,el)) ! calculate the normal of the interface in current ... + surfaceNormal = math_mul33x3(transpose(Fe(:,:,g,ip,el)), surfaceNormal_currentconf) / detFe ! ... and lattice configuration area = mesh_ipArea(n,ip,el) * math_norm3(surfaceNormal) + surfaceNormal = surfaceNormal / math_norm3(surfaceNormal) ! normalize the surface normal to unit length - ! normalize the surface normal to unit length - surfaceNormal = surfaceNormal / math_norm3(surfaceNormal) - - lineLength = 0.0_pReal - do s = 1,ns ! loop over slip systems - do t = 1,4 ! loop over dislocation types (only mobiles) - if ( sign(1.0_pReal,math_mul3x3(m(:,s,t),surfaceNormal)) == sign(1.0_pReal,gdot(s,t)) ) then - lineLength(s,t) = gdot(s,t) / constitutive_nonlocal_burgersPerSlipSystem(s,myInstance) & - * math_mul3x3(m(:,s,t), surfaceNormal) * area & - * constitutive_nonlocal_transmissivity(misorientation(4,n), misorientation(1:3,n)) ! dislocation line length that leaves this interface per second; weighted by a transmissivity factor - - thisRhoDotSgl(s,t) = thisRhoDotSgl(s,t) - lineLength(s,t) / mesh_ipVolume(ip,el) ! subtract dislocation density rate (= line length over volume) that leaves through an interface from my dotState ... - - if ( neighboring_el > 0 .and. neighboring_ip > 0 ) then - !$OMP CRITICAL (dotStateLock) - !$OMP FLUSH (dotState) - dotState(1,neighboring_ip,neighboring_el)%p((t-1)*ns+s) = dotState(1,neighboring_ip,neighboring_el)%p((t-1)*ns+s) & - + lineLength(s,t) / mesh_ipVolume(neighboring_ip,neighboring_el) ! ... and add it to the neighboring dotState (if neighbor exists) - !$OMP FLUSH (dotState) - !$OMPEND CRITICAL (dotStateLock) + do s = 1,ns + do t = 1,4 + if ( constitutive_nonlocal_v(s,t,g,ip,el) * math_mul3x3(m(:,s,t), surfaceNormal) > 0.0_pReal ) then ! outgoing flux + + if (neighboring_el > 0 .and. neighboring_ip > 0) then ! I have a neighbor ... + if (transmissivity > 0.99_pReal) then ! ... with perfect transmissivity of the interface + vAvg = 0.5_pReal * (constitutive_nonlocal_v(s,t,g,ip,el) + constitutive_nonlocal_v(s,t,g,neighboring_ip,neighboring_el)) + rhoAvg = 0.5_pReal * ( rhoSgl(s,t) + state(g,neighboring_ip,neighboring_el)%p((t-1)*ns+s) ) + else ! ... with low transmissivity of the interface, i.e. a discontinuity in the dislocation flux + vAvg = constitutive_nonlocal_v(s,t,g,ip,el) + rhoAvg = rhoSgl(s,t) + endif + + else ! no neighbor + vAvg = constitutive_nonlocal_v(s,t,g,ip,el) + rhoAvg = rhoSgl(s,t) endif + + lineLength = rhoAvg * vAvg * math_mul3x3(m(:,s,t), surfaceNormal) * area ! line length that wants to leave this interface + thisRhoDotSgl(s,t) = thisRhoDotSgl(s,t) - lineLength / mesh_ipVolume(ip,el) * transmissivity ! subtract positive dislocation flux that leaves the material point + thisRhoDotSgl(s,t+4) = thisRhoDotSgl(s,t+4) + lineLength / mesh_ipVolume(ip,el) * (1.0_pReal - transmissivity) & + * sign(1.0_pReal, vAvg) ! dislocation flux that is not able to leave through interface (because of low transmissivity) will remain as immobile single density at the material point + + else ! incoming flux + + if (neighboring_el > 0 .and. neighboring_ip > 0) then ! I have a neighbor ... + if (transmissivity > 0.99_pReal) then ! ... with perfect transmissivity of the interface + vAvg = 0.5_pReal * (constitutive_nonlocal_v(s,t,g,ip,el) + constitutive_nonlocal_v(s,t,g,neighboring_ip,neighboring_el)) + rhoAvg = 0.5_pReal * ( rhoSgl(s,t) + state(g,neighboring_ip,neighboring_el)%p((t-1)*ns+s) ) + else ! ... with low transmissivity of the interface, i.e. a discontinuity in the dislocation flux + vAvg = constitutive_nonlocal_v(s,t,g,neighboring_ip,neighboring_el) + rhoAvg = state(g,neighboring_ip,neighboring_el)%p((t-1)*ns+s) + endif + + else ! no neighbor + vAvg = 0.0_pReal + rhoAvg = 0.0_pReal + endif + + lineLength = rhoAvg * vAvg * math_mul3x3(m(:,s,t), surfaceNormal) * area ! line length that wants to leave this interface + thisRhoDotSgl(s,t) = thisRhoDotSgl(s,t) - lineLength / mesh_ipVolume(ip,el) * transmissivity ! subtract negative dislocation flux that enters the material point + endif enddo enddo + enddo totalRhoDotSgl = totalRhoDotSgl + thisRhoDotSgl totalRhoDotDip = totalRhoDotDip + thisRhoDotDip -if (debugger) then +if (selectiveDebugger) then !$OMP CRITICAL (write2out) - write(6,'(a,/,4(12(e12.5,x),/))') 'dislocation flux', thisRhoDotSgl(:,1:4) * timestep + write(6,'(a,/,8(12(e12.5,x),/))') 'dislocation flux', thisRhoDotSgl * timestep !$OMPEND CRITICAL (write2out) endif @@ -1401,18 +1484,32 @@ endif !*** formation by glide +thisRhoDotSgl = 0.0_pReal +thisRhoDotDip = 0.0_pReal + do c = 1,2 - thisRhoDotSgl(:,2*c-1) = - 4.0_pReal * dUpper(:,c) * rhoSgl(:,2*c-1) * rhoSgl(:,2*c) * abs(v) - thisRhoDotSgl(:,2*c) = thisRhoDotSgl(:,2*c-1) - thisRhoDotSgl(:,2*c+3) = - 2.0_pReal * dUpper(:,c) * abs(rhoSgl(:,2*c+3)) * rhoSgl(:,2*c) * abs(v) - thisRhoDotSgl(:,2*c+4) = - 2.0_pReal * dUpper(:,c) * abs(rhoSgl(:,2*c+4)) * rhoSgl(:,2*c-1) * abs(v) - thisRhoDotDip(:,c) = - thisRhoDotSgl(:,2*c-1) - thisRhoDotSgl(:,2*c) - thisRhoDotSgl(:,2*c+3) - thisRhoDotSgl(:,2*c+4) + + thisRhoDotSgl(:,2*c-1) = - 2.0_pReal * dUpper(:,c) / constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) & + * (rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) + rhoSgl(:,2*c) * abs(gdot(:,2*c-1)) & ! negative mobile <-> positive mobile + + abs(rhoSgl(:,2*c+4)) * abs(gdot(:,2*c-1))) ! negative immobile <-> positive mobile + + thisRhoDotSgl(:,2*c) = - 2.0_pReal * dUpper(:,c) / constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) & + * (rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) + rhoSgl(:,2*c) * abs(gdot(:,2*c-1)) & ! negative mobile <-> positive mobile + + abs(rhoSgl(:,2*c+3)) * abs(gdot(:,2*c))) ! negative mobile <-> positive immobile + + thisRhoDotSgl(:,2*c+3) = - 2.0_pReal * dUpper(:,c) / constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) & ! negative mobile <-> positive immobile + * rhoSgl(:,2*c+3) * abs(gdot(:,2*c)) + + thisRhoDotSgl(:,2*c+4) = - 2.0_pReal * dUpper(:,c) / constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) & + * rhoSgl(:,2*c+4) * abs(gdot(:,2*c-1)) ! negative immobile <-> positive mobile + + thisRhoDotDip(:,c) = - thisRhoDotSgl(:,2*c-1) - thisRhoDotSgl(:,2*c) + abs(thisRhoDotSgl(:,2*c+3)) + abs(thisRhoDotSgl(:,2*c+4)) enddo totalRhoDotSgl = totalRhoDotSgl + thisRhoDotSgl totalRhoDotDip = totalRhoDotDip + thisRhoDotDip -if (debugger) then +if (selectiveDebugger) then !$OMP CRITICAL (write2out) write(6,'(a,/,10(12(e12.5,x),/))') 'dipole formation by glide', thisRhoDotSgl * timestep, thisRhoDotDip * timestep !$OMPEND CRITICAL (write2out) @@ -1422,16 +1519,16 @@ endif !*** athermal annihilation forall (c=1:2) & - thisRhoDotDip(:,c) = - 2.0_pReal * dLower(:,c) * abs(v) & - * ( 4.0_pReal * rhoSgl(:,2*c-1) * rhoSgl(:,2*c) & ! was single hitting single - + abs(rhoSgl(:,2*c+3)) * rhoSgl(:,2*c) + abs(rhoSgl(:,2*c+4)) * rhoSgl(:,2*c-1) & ! was single hitting immobile/used single - + rhoDip(:,c) * ( rhoSgl(:,2*c-1) + rhoSgl(:,2*c) ) ) ! single knocks dipole constituent -thisRhoDotSgl = 0.0_pReal ! singles themselves don't annihilate + thisRhoDotDip(:,c) = - 2.0_pReal * dLower(:,c) / constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) & + * ( 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/used single + + rhoDip(:,c) * ( abs(gdot(:,2*c-1)) + abs(gdot(:,2*c)) ) ) ! single knocks dipole constituent +thisRhoDotSgl = 0.0_pReal ! singles themselves don't annihilate totalRhoDotSgl = totalRhoDotSgl + thisRhoDotSgl totalRhoDotDip = totalRhoDotDip + thisRhoDotDip -if (debugger) then +if (selectiveDebugger) then !$OMP CRITICAL (write2out) write(6,'(a,/,2(12(e12.5,x),/))') 'athermal dipole annihilation', thisRhoDotDip * timestep !$OMPEND CRITICAL (write2out) @@ -1446,14 +1543,14 @@ vClimb = constitutive_nonlocal_atomicVolume(myInstance) * D / ( kB * Temperatur * constitutive_nonlocal_Gmod(myInstance) / ( 2.0_pReal * pi * (1.0_pReal-constitutive_nonlocal_nu(myInstance)) ) & * 2.0_pReal / ( dUpper(:,1) + dLower(:,1) ) -thisRhoDotDip(:,1) = - 4.0_pReal * rhoDip(:,1) * vClimb / ( dUpper(:,1) - dLower(:,1) ) ! edge climb -thisRhoDotDip(:,2) = 0.0_pReal !!! cross slipping still has to be implemented !!! +thisRhoDotDip(:,1) = - 4.0_pReal * rhoDip(:,1) * vClimb / ( dUpper(:,1) - dLower(:,1) ) ! edge climb +thisRhoDotDip(:,2) = 0.0_pReal !!! cross slipping still has to be implemented !!! thisRhoDotSgl = 0.0_pReal totalRhoDotSgl = totalRhoDotSgl + thisRhoDotSgl totalRhoDotDip = totalRhoDotDip + thisRhoDotDip -if (debugger) then +if (selectiveDebugger) then !$OMP CRITICAL (write2out) write(6,'(a,/,2(12(e12.5,x),/))') 'thermally activated dipole annihilation', thisRhoDotDip * timestep !$OMPEND CRITICAL (write2out) @@ -1489,7 +1586,7 @@ endif dotState(1,ip,el)%p(1:8*ns) = dotState(1,ip,el)%p(1:8*ns) + reshape(totalRhoDotSgl,(/8*ns/)) dotState(1,ip,el)%p(8*ns+1:10*ns) = dotState(1,ip,el)%p(8*ns+1:10*ns) + reshape(totalRhoDotDip,(/2*ns/)) -if (debugger) then +if (selectiveDebugger) then !$OMP CRITICAL (write2out) write(6,'(a,/,8(12(e12.5,x),/))') 'deltaRho:', totalRhoDotSgl * timestep write(6,'(a,/,2(12(e12.5,x),/))') 'deltaRhoDip:', totalRhoDotDip * timestep @@ -1653,10 +1750,9 @@ real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstan lineLength ! dislocation line length leaving the current interface real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: & rhoForest, & ! forest dislocation density - tauSlipThreshold, & ! threshold shear stress - tauSlip, & ! current resolved shear stress - previousTauSlip, & ! previous resolved shear stress - v, & ! dislocation velocity + tauThreshold, & ! threshold shear stress + tau, & ! current resolved shear stress + previousTau, & ! previous resolved shear stress invLambda, & ! inverse of mean free path for dislocations vClimb ! climb velocity of edge dipoles real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),2) :: & @@ -1667,8 +1763,9 @@ real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstan dUpper, & ! current maximum stable dipole distance for edges and screws previousDUpper, & ! previous maximum stable dipole distance for edges and screws dUpperDot ! rate of change of the maximum stable dipole distance for edges and screws -real(pReal), dimension(3,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4) :: & - m ! direction of dislocation motion +real(pReal), dimension(3,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),2) :: & + m, & ! direction of dislocation motion for edge and screw (unit vector) + m_currentconf ! direction of dislocation motion for edge and screw (unit vector) in current configuration real(pReal), dimension(3,3) :: F, & ! total deformation gradient neighboring_F, & ! total deformation gradient of my neighbor Favg ! average total deformation gradient of me and my neighbor @@ -1696,7 +1793,7 @@ forall (t = 1:8) previousRhoSgl(:,t) = previousState(g,ip,el)%p((t-1)*ns+1:t*ns) forall (c = 1:2) rhoDip(:,c) = state(g,ip,el)%p((7+c)*ns+1:(8+c)*ns) forall (c = 1:2) previousRhoDip(:,c) = previousState(g,ip,el)%p((7+c)*ns+1:(8+c)*ns) rhoForest = state(g,ip,el)%p(10*ns+1:11*ns) -tauSlipThreshold = state(g,ip,el)%p(11*ns+1:12*ns) +tauThreshold = state(g,ip,el)%p(11*ns+1:12*ns) Tdislocation_v = state(g,ip,el)%p(12*ns+1:12*ns+6) previousTdislocation_v = previousState(g,ip,el)%p(12*ns+1:12*ns+6) forall (t = 1:8) rhoDotSgl(:,t) = dotState(g,ip,el)%p((t-1)*ns+1:t*ns) @@ -1705,38 +1802,34 @@ forall (c = 1:2) rhoDotDip(:,c) = dotState(g,ip,el)%p((7+c)*ns+1:(8+c)*ns) !* Calculate shear rate -do s = 1,ns - sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance) - tauSlip(s) = math_mul6x6( Tstar_v + Tdislocation_v, lattice_Sslip_v(:,sLattice,myStructure) ) - previousTauSlip(s) = math_mul6x6( previousTstar_v + previousTdislocation_v, lattice_Sslip_v(:,sLattice,myStructure) ) -enddo - -v = constitutive_nonlocal_v0PerSlipSystem(:,myInstance) & - * exp( - constitutive_nonlocal_Q0(myInstance) / ( kB * Temperature) * (1.0_pReal - (abs(tauSlip)/tauSlipThreshold) ) ) & - * sign(1.0_pReal,tauSlip) - do t = 1,4 do s = 1,ns - if (rhoSgl(s,t+4) * v(s) < 0.0_pReal) then + if (rhoSgl(s,t+4) * constitutive_nonlocal_v(s,t,g,ip,el) < 0.0_pReal) then rhoSgl(s,t) = rhoSgl(s,t) + abs(rhoSgl(s,t+4)) ! remobilization of immobile singles for changing sign of v (bauschinger effect) rhoSgl(s,t+4) = 0.0_pReal ! remobilization of immobile singles for changing sign of v (bauschinger effect) endif enddo enddo -forall (t = 1:4) gdot(:,t) = rhoSgl(:,t) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) * v - +forall (t = 1:4) & + gdot(:,t) = rhoSgl(:,t) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) * constitutive_nonlocal_v(:,t,g,ip,el) !* calculate limits for stable dipole height and its rate of change +do s = 1,ns + sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance) + tau(s) = math_mul6x6( Tstar_v + Tdislocation_v, lattice_Sslip_v(:,sLattice,myStructure) ) + previousTau(s) = math_mul6x6( previousTstar_v + previousTdislocation_v, lattice_Sslip_v(:,sLattice,myStructure) ) +enddo + dLower(:,1) = constitutive_nonlocal_dLowerEdgePerSlipSystem(:,myInstance) dLower(:,2) = constitutive_nonlocal_dLowerScrewPerSlipSystem(:,myInstance) dUpper(:,2) = min( constitutive_nonlocal_Gmod(myInstance) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) & - / ( 8.0_pReal * pi * abs(tauSlip) ), & + / ( 8.0_pReal * pi * abs(tau) ), & 1.0_pReal / sqrt( sum(abs(rhoSgl),2)+sum(rhoDip,2) ) ) dUpper(:,1) = dUpper(:,2) / ( 1.0_pReal - constitutive_nonlocal_nu(myInstance) ) previousDUpper(:,2) = min( constitutive_nonlocal_Gmod(myInstance) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) & - / ( 8.0_pReal * pi * abs(previousTauSlip) ), & + / ( 8.0_pReal * pi * abs(previousTau) ), & 1.0_pReal / sqrt( sum(abs(previousRhoSgl),2) + sum(previousRhoDip,2) ) ) previousDUpper(:,1) = previousDUpper(:,2) / ( 1.0_pReal - constitutive_nonlocal_nu(myInstance) ) @@ -1747,47 +1840,12 @@ else endif -!* calculate fluxes - -m(:,:,1) = lattice_sd(:, constitutive_nonlocal_slipSystemLattice(:,myInstance), myStructure) -m(:,:,2) = -lattice_sd(:, constitutive_nonlocal_slipSystemLattice(:,myInstance), myStructure) -m(:,:,3) = lattice_st(:, constitutive_nonlocal_slipSystemLattice(:,myInstance), myStructure) -m(:,:,4) = -lattice_st(:, constitutive_nonlocal_slipSystemLattice(:,myInstance), myStructure) - -F = math_mul33x33(Fe(:,:,g,ip,el), Fp(:,:,g,ip,el)) -detFe = math_det3x3(Fe(:,:,g,ip,el)) - -fluxes = 0.0_pReal -do n = 1,FE_NipNeighbors(mesh_element(2,el)) ! loop through my neighbors - - neighboring_el = mesh_ipNeighborhood(1,n,ip,el) - neighboring_ip = mesh_ipNeighborhood(2,n,ip,el) - - if ( neighboring_el > 0 .and. neighboring_ip > 0 ) then ! if neighbor exists, total deformation gradient is averaged over me and my neighbor - neighboring_F = math_mul33x33(Fe(:,:,g,neighboring_ip,neighboring_el), Fp(:,:,g,neighboring_ip,neighboring_el)) - Favg = 0.5_pReal * (F + neighboring_F) - else - Favg = F - endif - - surfaceNormal_currentconf = math_det3x3(Favg) * math_mul33x3(math_inv3x3(transpose(Favg)), mesh_ipAreaNormal(:,n,ip,el)) ! normal of the interface in current configuration - surfaceNormal = math_mul33x3(transpose(Fe(:,:,g,ip,el)), surfaceNormal_currentconf) / detFe ! normal of the interface in lattice configuration - surfaceNormal = surfaceNormal / math_norm3(surfaceNormal) ! normalize the surface normal to unit length - - area = mesh_ipArea(n,ip,el) * math_norm3(surfaceNormal) ! area of the interface - - do s = 1,ns ! loop over slip systems - do t = 1,4 ! loop over dislocation types - if ( sign(1.0_pReal,math_mul3x3(m(:,s,t),surfaceNormal)) == sign(1.0_pReal,gdot(s,t)) ) & - fluxes(s,n,t) = gdot(s,t) / constitutive_nonlocal_burgersPerSlipSystem(s,myInstance) & - * math_mul3x3(m(:,s,t), surfaceNormal) * area & - * constitutive_nonlocal_transmissivity(misorientation(4,n), misorientation(1:3,n)) & - / mesh_ipVolume(ip,el) - enddo - enddo -enddo - +!*** dislocation motion +m(:,:,1) = lattice_sd(:, constitutive_nonlocal_slipSystemLattice(:,myInstance), myStructure) +m(:,:,2) = lattice_st(:, constitutive_nonlocal_slipSystemLattice(:,myInstance), myStructure) +forall (c = 1:2, s = 1:ns) & + m_currentconf(:,s,c) = math_mul33x3(Fe(:,:,g,ip,el), m(:,s,c)) do o = 1,phase_Noutput(material_phase(g,ip,el)) @@ -1941,8 +1999,22 @@ do o = 1,phase_Noutput(material_phase(g,ip,el)) enddo cs = cs + ns + case ('resolvedstress_internal') + do s = 1,ns + sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance) + constitutive_nonlocal_postResults(cs+s) = math_mul6x6(Tdislocation_v, lattice_Sslip_v(:,sLattice,myStructure) ) + enddo + cs = cs + ns + + case ('resolvedstress_external') + do s = 1,ns + sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance) + constitutive_nonlocal_postResults(cs+s) = math_mul6x6(Tstar_v, lattice_Sslip_v(:,sLattice,myStructure) ) + enddo + cs = cs + ns + case ('resistance') - constitutive_nonlocal_postResults(cs+1:cs+ns) = tauSlipThreshold + constitutive_nonlocal_postResults(cs+1:cs+ns) = tauThreshold cs = cs + ns case ('rho_dot') @@ -1978,8 +2050,9 @@ do o = 1,phase_Noutput(material_phase(g,ip,el)) case ('rho_dot_sgl2dip') do c=1,2 ! dipole formation by glide constitutive_nonlocal_postResults(cs+1:cs+ns) = constitutive_nonlocal_postResults(cs+1:cs+ns) + & - 2.0_pReal * dUpper(:,c) * abs(v) * ( 4.0_pReal * rhoSgl(:,2*c-1) * rhoSgl(:,2*c) & ! was single hitting single - + abs(rhoSgl(:,2*c+3)) * rhoSgl(:,2*c) + abs(rhoSgl(:,2*c+4)) * rhoSgl(:,2*c-1) ) ! was single hitting immobile/used single + 2.0_pReal * dUpper(:,c) / constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) & + * ( 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/used single enddo ! do c=1,2 ! forall (s=1:ns, dUpperDot(s,c) > 0.0_pReal) & ! dipole formation by stress decrease @@ -2000,9 +2073,10 @@ do o = 1,phase_Noutput(material_phase(g,ip,el)) case ('rho_dot_ann_ath') do c=1,2 constitutive_nonlocal_postResults(cs+1:cs+ns) = constitutive_nonlocal_postResults(cs+1:cs+ns) + & - 2.0_pReal * dLower(:,c) * abs(v) * ( 4.0_pReal * rhoSgl(:,2*c-1) * rhoSgl(:,2*c) & ! was single hitting single - + abs(rhoSgl(:,2*c+3)) * rhoSgl(:,2*c) + abs(rhoSgl(:,2*c+4)) * rhoSgl(:,2*c-1) & ! was single hitting immobile/used single - + rhoDip(:,c) * ( rhoSgl(:,2*c-1) + rhoSgl(:,2*c) ) ) ! single knocks dipole constituent + 2.0_pReal * dLower(:,c) / constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) & + * ( 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/used single + + rhoDip(:,c) * ( rhoSgl(:,2*c-1) + rhoSgl(:,2*c) ) ) ! single knocks dipole constituent enddo cs = cs + ns @@ -2017,121 +2091,52 @@ do o = 1,phase_Noutput(material_phase(g,ip,el)) ! !!! cross-slip of screws missing !!! cs = cs + ns - case ('rho_dot_flux_n1_ep') - constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,1,1) + case ('fluxdensity_edge_pos_x') + constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(:,1) * constitutive_nonlocal_v(:,1,g,ip,el) * m_currentconf(1,:,1) cs = cs + ns - case ('rho_dot_flux_n1_en') - constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,1,2) + case ('fluxdensity_edge_pos_y') + constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(:,1) * constitutive_nonlocal_v(:,1,g,ip,el) * m_currentconf(2,:,1) cs = cs + ns - case ('rho_dot_flux_n1_sp') - constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,1,3) + case ('fluxdensity_edge_pos_z') + constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(:,1) * constitutive_nonlocal_v(:,1,g,ip,el) * m_currentconf(3,:,1) cs = cs + ns - case ('rho_dot_flux_n1_sn') - constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,1,4) + case ('fluxdensity_edge_neg_x') + constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(:,2) * constitutive_nonlocal_v(:,2,g,ip,el) * m_currentconf(1,:,2) cs = cs + ns - case ('rho_dot_flux_n2_ep') - constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,2,1) + case ('fluxdensity_edge_neg_y') + constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(:,2) * constitutive_nonlocal_v(:,2,g,ip,el) * m_currentconf(2,:,2) cs = cs + ns - case ('rho_dot_flux_n2_en') - constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,2,2) + case ('fluxdensity_edge_neg_z') + constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(:,2) * constitutive_nonlocal_v(:,2,g,ip,el) * m_currentconf(3,:,2) cs = cs + ns - case ('rho_dot_flux_n2_sp') - constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,2,3) + case ('fluxdensity_screw_pos_x') + constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(:,3) * constitutive_nonlocal_v(:,3,g,ip,el) * m_currentconf(1,:,3) cs = cs + ns - case ('rho_dot_flux_n2_sn') - constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,2,4) + case ('fluxdensity_screw_pos_y') + constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(:,3) * constitutive_nonlocal_v(:,3,g,ip,el) * m_currentconf(2,:,3) cs = cs + ns - case ('rho_dot_flux_n3_ep') - constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,3,1) + case ('fluxdensity_screw_pos_z') + constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(:,3) * constitutive_nonlocal_v(:,3,g,ip,el) * m_currentconf(3,:,3) cs = cs + ns - case ('rho_dot_flux_n3_en') - constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,3,2) + case ('fluxdensity_screw_neg_x') + constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(:,4) * constitutive_nonlocal_v(:,4,g,ip,el) * m_currentconf(1,:,4) cs = cs + ns - case ('rho_dot_flux_n3_sp') - constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,3,3) + case ('fluxdensity_screw_neg_y') + constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(:,4) * constitutive_nonlocal_v(:,4,g,ip,el) * m_currentconf(2,:,4) cs = cs + ns - case ('rho_dot_flux_n3_sn') - constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,3,4) - cs = cs + ns - - case ('rho_dot_flux_n4_ep') - constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,4,1) - cs = cs + ns - - case ('rho_dot_flux_n4_en') - constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,4,2) - cs = cs + ns - - case ('rho_dot_flux_n4_sp') - constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,4,3) - cs = cs + ns - - case ('rho_dot_flux_n4_sn') - constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,4,4) - cs = cs + ns - - case ('rho_dot_flux_n5_ep') - constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,5,1) - cs = cs + ns - - case ('rho_dot_flux_n5_en') - constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,5,2) - cs = cs + ns - - case ('rho_dot_flux_n5_sp') - constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,5,3) - cs = cs + ns - - case ('rho_dot_flux_n5_sn') - constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,5,4) - cs = cs + ns - - case ('rho_dot_flux_n6_ep') - constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,6,1) - cs = cs + ns - - case ('rho_dot_flux_n6_en') - constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,6,2) - cs = cs + ns - - case ('rho_dot_flux_n6_sp') - constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,6,3) - cs = cs + ns - - case ('rho_dot_flux_n6_sn') - constitutive_nonlocal_postResults(cs+1:cs+ns) = fluxes(:,6,4) - cs = cs + ns - - case ('rho_dot_flux_ep') - constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(fluxes(:,:,1), 2) - cs = cs + ns - - case ('rho_dot_flux_en') - constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(fluxes(:,:,2), 2) - cs = cs + ns - - case ('rho_dot_flux_sp') - constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(fluxes(:,:,3), 2) - cs = cs + ns - - case ('rho_dot_flux_sn') - constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(fluxes(:,:,4), 2) - cs = cs + ns - - case ('rho_dot_flux') - constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(fluxes(:,:,1), 2) + sum(fluxes(:,:,2), 2) & - + sum(fluxes(:,:,3), 2) + sum(fluxes(:,:,4), 2) + case ('fluxdensity_screw_neg_z') + constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(:,4) * constitutive_nonlocal_v(:,4,g,ip,el) * m_currentconf(3,:,4) cs = cs + ns case ('d_upper_edge') diff --git a/code/crystallite.f90 b/code/crystallite.f90 index c0f72b52e..629d51b39 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -249,6 +249,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco) nState, & nCryst use debug, only: debugger, & + selectiveDebugger, & + debug_e, & + debug_i, & + debug_g, & debug_CrystalliteLoopDistribution, & debug_CrystalliteStateLoopDistribution, & debug_StiffnessStateLoopDistribution @@ -338,6 +342,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) logical, dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & storedConvergenceFlag logical, dimension(3,3) :: mask + logical forceLocalStiffnessCalculation = .true. ! flag indicating that stiffness calculation is always done locally ! ------ initialize to starting condition ------ @@ -388,9 +393,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco) myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains -! debugger = (e == 1 .and. i == 1 .and. g == 1) + selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) if (crystallite_converged(g,i,e)) then - if (debugger) then + if (selectiveDebugger) then !$OMP CRITICAL (write2out) write(6,'(a21,f10.8,a32,f10.8,a35)') 'winding forward from ', & crystallite_subFrac(g,i,e),' to current crystallite_subfrac ', & @@ -422,7 +427,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) crystallite_Lp(:,:,g,i,e) = crystallite_subLp0(:,:,g,i,e) ! ...plastic velocity grad constitutive_state(g,i,e)%p = constitutive_subState0(g,i,e)%p ! ...microstructure crystallite_Tstar_v(:,g,i,e) = crystallite_subTstar0_v(:,g,i,e) ! ...2nd PK stress - if (debugger) then + if (selectiveDebugger) then !$OMP CRITICAL (write2out) write(6,'(a78,f10.8)') 'cutback step in crystallite_stressAndItsTangent with new crystallite_subStep: ',& crystallite_subStep(g,i,e) @@ -463,7 +468,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains if (crystallite_todo(g,i,e)) then ! all undone crystallites - call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states + call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, & + crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states constitutive_previousDotState2(g,i,e)%p = 0.0_pReal constitutive_previousDotState(g,i,e)%p = 0.0_pReal constitutive_dotState(g,i,e)%p = 0.0_pReal ! zero out dotStates @@ -475,7 +481,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains - ! debugger = (e == 1 .and. i == 1 .and. g == 1) +! selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) if (crystallite_todo(g,i,e)) then ! all undone crystallites call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & @@ -491,7 +497,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains - ! debugger = (e == 1 .and. i == 1 .and. g == 1) +! selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) if (crystallite_todo(g,i,e)) then ! all undone crystallites crystallite_stateConverged(g,i,e) = crystallite_updateState(g,i,e) ! update state crystallite_temperatureConverged(g,i,e) = crystallite_updateTemperature(g,i,e) ! update temperature @@ -520,7 +526,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains - ! debugger = (e == 1 .and. i == 1 .and. g == 1) +! selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) if (crystallite_todo(g,i,e)) & ! all undone crystallites crystallite_onTrack(g,i,e) = crystallite_integrateStress(g,i,e) enddo; enddo; enddo @@ -571,7 +577,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains - ! debugger = (e == 1 .and. i == 1 .and. g == 1) +! selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) if (crystallite_todo(g,i,e)) then ! all undone crystallites call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & @@ -594,7 +600,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains - ! debugger = (e == 1 .and. i == 1 .and. g == 1) +! selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) if (crystallite_todo(g,i,e)) then ! all undone crystallites crystallite_stateConverged(g,i,e) = crystallite_updateState(g,i,e) ! update state crystallite_temperatureConverged(g,i,e) = crystallite_updateTemperature(g,i,e) ! update temperature @@ -615,7 +621,11 @@ subroutine crystallite_stressAndItsTangent(updateJaco) !$OMP CRITICAL (write2out) write(6,*) count(crystallite_converged(:,:,:)),'grains converged after state integration no.', NiterationState write(6,*) - !write(6,'(8(L,x))') crystallite_converged(:,:,:) +! write(6,'(8(L,x))') crystallite_converged(:,:,:) +! do e = FEsolving_execElem(1),FEsolving_execElem(2) +! if (any(.not. crystallite_converged(:,:,e))) & +! write(6,'(i4,8(x,L))') e, crystallite_converged(:,:,e) +! enddo !$OMPEND CRITICAL (write2out) endif if (any(.not. crystallite_converged .and. .not. crystallite_localConstitution)) & ! any non-local not yet converged? @@ -684,24 +694,24 @@ subroutine crystallite_stressAndItsTangent(updateJaco) storedP = crystallite_P storedConvergenceFlag = crystallite_converged - if (all(crystallite_localConstitution) .or. theInc < 2) then ! all grains have local constitution, so local convergence of perturbed grain is sufficient + if (all(crystallite_localConstitution) .or. theInc < 2 .or. forceLocalStiffnessCalculation) then ! all grains have local constitution, so local convergence of perturbed grain is sufficient !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains -! debugger = (e==1 .and. i==1 .and. g==1) +! selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) if (crystallite_requested(g,i,e)) then ! first check whether is requested at all! if (crystallite_converged(g,i,e)) then ! grain converged in above iteration - if (debugger) then + if (selectiveDebugger) then !$OMP CRITICAL (write2out) write (6,*) '#############' write (6,*) 'central solution of cryst_StressAndTangent' write (6,*) '#############' - write (6,'(a,/,3(3(f12.4,x)/))') ' P of 1 1 1',storedP(1:3,:,g,i,e)/1e6 - write (6,'(a,/,3(3(f12.8,x)/))') ' Fp of 1 1 1',storedFp(1:3,:,g,i,e) - write (6,'(a,/,3(3(f12.8,x)/))') ' Lp of 1 1 1',storedLp(1:3,:,g,i,e) + write (6,'(a8,3(x,i4),/,3(3(f12.4,x)/))') ' P of', g, i, e, storedP(1:3,:,g,i,e)/1e6 + write (6,'(a8,3(x,i4),/,3(3(f12.8,x)/))') ' Fp of', g, i, e, storedFp(1:3,:,g,i,e) + write (6,'(a8,3(x,i4),/,3(3(f12.8,x)/))') ' Lp of', g, i, e, storedLp(1:3,:,g,i,e) !$OMPEND CRITICAL (write2out) endif @@ -711,10 +721,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco) do k = 1,3 ! perturbation... do l = 1,3 ! ...components to the positive direction crystallite_subF(k,l,g,i,e) = crystallite_subF(k,l,g,i,e) + myPert ! perturb single component (either forward or backward) - if (debugger) then + if (selectiveDebugger) then !$OMP CRITICAL (write2out) write (6,'(i1,x,i1)') k,l - write (6,'(a,/,3(3(f12.6,x)/))') 'pertF of 1 1 1',crystallite_subF(1:3,:,g,i,e) + write (6,'(a8,3(x,i4),/,3(3(f12.6,x)/))') 'pertF of', g, i, e, crystallite_subF(1:3,:,g,i,e) !$OMPEND CRITICAL (write2out) endif onTrack = .true. @@ -733,12 +743,13 @@ subroutine crystallite_stressAndItsTangent(updateJaco) temperatureConverged = crystallite_updateTemperature(g,i,e) ! update temperature converged = stateConverged .and. temperatureConverged endif - if (debugger) then + if (selectiveDebugger) then !$OMP CRITICAL (write2out) write (6,*) '-------------' write (6,'(a,x,l,x,l)') 'ontrack + converged:',onTrack,converged - write (6,'(a,/,3(3(f12.4,x)/))') 'pertP/MPa of 1 1 1',crystallite_P(1:3,:,g,i,e)/1e6 - write (6,'(a,/,3(3(f12.4,x)/))') 'DP/MPa of 1 1 1',(crystallite_P(1:3,:,g,i,e)-storedP(1:3,:,g,i,e))/1e6 + write (6,'(a12,3(x,i4),/,3(3(f12.4,x)/))') 'pertP/MPa of', g, i, e, crystallite_P(1:3,:,g,i,e)/1e6 + write (6,'(a12,3(x,i4),/,3(3(f12.4,x)/))') 'DP/MPa of', g, i, e, & + (crystallite_P(1:3,:,g,i,e)-storedP(1:3,:,g,i,e))/1e6 !$OMPEND CRITICAL (write2out) endif enddo @@ -918,7 +929,8 @@ endsubroutine constitutive_state, & constitutive_relevantState, & constitutive_microstructure - use debug, only: debugger + use debug, only: debugger, & + selectiveDebugger use FEsolving, only: cycleCounter, theInc !*** input variables ***! @@ -957,13 +969,14 @@ endsubroutine ! update the microstructure constitutive_state(g,i,e)%p(1:mySize) = constitutive_state(g,i,e)%p(1:mySize) - residuum - call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe, crystallite_Fp, g, i, e) + call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, crystallite_Fp, & + g, i, e) ! setting flag to true if state is below relative tolerance, otherwise set it to false crystallite_updateState = all( constitutive_state(g,i,e)%p(1:mySize) < constitutive_relevantState(g,i,e)%p(1:mySize) & .or. abs(residuum) < rTol_crystalliteState*abs(constitutive_state(g,i,e)%p(1:mySize))) - if (debugger) then + if (selectiveDebugger) then !$OMP CRITICAL (write2out) if (crystallite_updateState) then write(6,*) '::: updateState converged',g,i,e @@ -1073,6 +1086,7 @@ endsubroutine iJacoLpresiduum, & relevantStrain use debug, only: debugger, & + selectiveDebugger, & debug_cumLpCalls, & debug_cumLpTicks, & debug_StressLoopDistribution @@ -1214,13 +1228,12 @@ LpLoop: do debug_cumLpCalls = debug_cumLpCalls + 1_pInt debug_cumLpTicks = debug_cumLpTicks + tock-tick if (tock < tick) debug_cumLpTicks = debug_cumLpTicks + maxticks - if (debugger) then + if (selectiveDebugger) then !$OMP CRITICAL (write2out) write(6,*) '::: integrateStress at ' ,g,i,e, ' ; iteration ', NiterationStress write(6,*) write(6,'(a,/,3(3(f20.7,x)/))') 'Lp_constitutive', Lp_constitutive write(6,'(a,/,3(3(f20.7,x)/))') 'Lpguess', Lpguess - ! call flush(6) !$OMPEND CRITICAL (write2out) endif @@ -1332,7 +1345,7 @@ LpLoop: do ! set return flag to true crystallite_integrateStress = .true. - if (debugger) then + if (selectiveDebugger) then !$OMP CRITICAL (write2out) write(6,*) '::: integrateStress converged at ',g,i,e,' ; iteration ', NiterationStress write(6,*) @@ -1468,11 +1481,12 @@ logical error crystallite_R(:,:,1,neighboring_i,neighboring_e), & symmetryType) ! calculate misorientation - else - crystallite_misorientation(4,n,1,i,e) = 400.0_pReal -! call IO_warning(-1, e, i, g, 'couldnt calculate misorientation because of two different lattice structures') + else ! for neighbor with different phase + crystallite_misorientation(4,n,1,i,e) = 400.0_pReal ! set misorientation angle to 400 endif + else ! no existing neighbor + crystallite_misorientation(4,n,1,i,e) = 0.0_pReal ! set misorientation angle to zero endif enddo endif diff --git a/code/debug.f90 b/code/debug.f90 index c2e2acf3d..d2b63c7ed 100644 --- a/code/debug.f90 +++ b/code/debug.f90 @@ -17,6 +17,10 @@ integer(pInt) :: debug_cumLpCalls = 0_pInt integer(pInt) :: debug_cumDotStateCalls = 0_pInt integer(pInt) :: debug_cumDotTemperatureCalls = 0_pInt + integer(pInt) :: debug_e = 1_pInt + integer(pInt) :: debug_i = 1_pInt + integer(pInt) :: debug_g = 1_pInt + logical :: selectiveDebugger = .false. logical :: debugger = .false. logical :: distribution_init = .false. diff --git a/code/homogenization.f90 b/code/homogenization.f90 index cc05d8299..7f9d4c0ba 100644 --- a/code/homogenization.f90 +++ b/code/homogenization.f90 @@ -220,6 +220,9 @@ subroutine materialpoint_stressAndItsTangent(& crystallite_stressAndItsTangent, & crystallite_orientations use debug, only: debugger, & + selectiveDebugger, & + debug_e, & + debug_i, & debug_MaterialpointLoopDistribution, & debug_MaterialpointStateLoopDistribution @@ -278,11 +281,11 @@ subroutine materialpoint_stressAndItsTangent(& myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed - debugger = (e == 1 .and. i == 1) + selectiveDebugger = (e == debug_e .and. i == debug_i) ! if our materialpoint converged or consists of only one single grain then we are either finished or have to wind forward if ( materialpoint_converged(i,e) .or. (myNgrains == 1_pInt .and. materialpoint_subStep(i,e) <= 1.0_pReal) ) then - if (debugger) then + if (selectiveDebugger) then !$OMP CRITICAL (write2out) write(6,'(a21,f10.8,a34,f10.8,a37,/)') 'winding forward from ', & materialpoint_subFrac(i,e), ' to current materialpoint_subFrac ', & @@ -321,7 +324,7 @@ subroutine materialpoint_stressAndItsTangent(& materialpoint_subStep(i,e) = subStepSizeHomog * materialpoint_subStep(i,e) ! crystallite had severe trouble, so do a significant cutback ! <> - if (debugger) then + if (selectiveDebugger) then !$OMP CRITICAL (write2out) write(6,'(a82,f10.8,/)') 'cutback step in materialpoint_stressAndItsTangent with new materialpoint_subStep: ',& materialpoint_subStep(i,e) diff --git a/code/material.config b/code/material.config index fc9a83af1..935fc84de 100644 --- a/code/material.config +++ b/code/material.config @@ -149,6 +149,8 @@ constitution nonlocal (output) delta_dip (output) shearrate (output) resolvedstress +(output) resolvedstress_internal +(output) resolvedstress_external (output) resistance (output) rho_dot (output) rho_dot_sgl @@ -160,35 +162,18 @@ constitution nonlocal (output) rho_dot_dip2sgl (output) rho_dot_ann_ath (output) rho_dot_ann_the -(output) rho_dot_flux -(output) rho_dot_flux_ep -(output) rho_dot_flux_en -(output) rho_dot_flux_sp -(output) rho_dot_flux_sn -(output) rho_dot_flux_n1_ep -(output) rho_dot_flux_n1_en -(output) rho_dot_flux_n1_sp -(output) rho_dot_flux_n1_sn -(output) rho_dot_flux_n2_ep -(output) rho_dot_flux_n2_en -(output) rho_dot_flux_n2_sp -(output) rho_dot_flux_n2_sn -(output) rho_dot_flux_n3_ep -(output) rho_dot_flux_n3_en -(output) rho_dot_flux_n3_sp -(output) rho_dot_flux_n3_sn -(output) rho_dot_flux_n4_ep -(output) rho_dot_flux_n4_en -(output) rho_dot_flux_n4_sp -(output) rho_dot_flux_n4_sn -(output) rho_dot_flux_n5_ep -(output) rho_dot_flux_n5_en -(output) rho_dot_flux_n5_sp -(output) rho_dot_flux_n5_sn -(output) rho_dot_flux_n6_ep -(output) rho_dot_flux_n6_en -(output) rho_dot_flux_n6_sp -(output) rho_dot_flux_n6_sn +(output) fluxDensity_edge_pos_x +(output) fluxDensity_edge_pos_y +(output) fluxDensity_edge_pos_z +(output) fluxDensity_edge_neg_x +(output) fluxDensity_edge_neg_y +(output) fluxDensity_edge_neg_z +(output) fluxDensity_screw_pos_x +(output) fluxDensity_screw_pos_y +(output) fluxDensity_screw_pos_z +(output) fluxDensity_screw_neg_x +(output) fluxDensity_screw_neg_y +(output) fluxDensity_screw_neg_z (output) d_upper_edge (output) d_upper_screw (output) d_upper_dot_edge