From c7e3ac28f63bfae9a8aeb8f9232a975793c847cd Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 7 Feb 2020 11:44:03 +0100 Subject: [PATCH 1/7] preparing data handling for explicit forward of flux --- src/constitutive_plastic_nonlocal.f90 | 615 ++++++++++++++------------ 1 file changed, 321 insertions(+), 294 deletions(-) diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index 42b8a34f5..358fd964d 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -11,7 +11,7 @@ submodule(constitutive) plastic_nonlocal IPvolume => geometry_plastic_nonlocal_IPvolume0, & IParea => geometry_plastic_nonlocal_IParea0, & IPareaNormal => geometry_plastic_nonlocal_IPareaNormal0 - + real(pReal), parameter :: & KB = 1.38e-23_pReal !< Physical parameter, Boltzmann constant in J/Kelvin @@ -44,11 +44,11 @@ submodule(constitutive) plastic_nonlocal integer, dimension(:), allocatable :: & totalNslip !< total number of active slip systems for each instance !END DEPRECATED - + real(pReal), dimension(:,:,:,:,:,:), allocatable :: & compatibility !< slip system compatibility between me and my neighbors - - enum, bind(c) + + enum, bind(c) enumerator :: & undefined_ID, & rho_sgl_mob_edg_pos_ID, & @@ -73,7 +73,7 @@ submodule(constitutive) plastic_nonlocal v_scr_neg_ID, & gamma_ID end enum - + type :: tParameters !< container type for internal constitutive parameters real(pReal) :: & atomicVolume, & !< atomic volume @@ -135,27 +135,27 @@ submodule(constitutive) plastic_nonlocal integer, dimension(:) ,allocatable :: & Nslip,& colinearSystem !< colinear system to the active slip system (only valid for fcc!) - + logical :: & shortRangeStressCorrection, & !< flag indicating the use of the short range stress correction by a excess density gradient term probabilisticMultiplication - + integer(kind(undefined_ID)), dimension(:), allocatable :: & outputID !< ID of each post result output - + end type tParameters - + type :: tNonlocalMicrostructure real(pReal), allocatable, dimension(:,:) :: & - tau_pass, & - tau_Back + tau_pass, & + tau_Back end type tNonlocalMicrostructure - + type :: tNonlocalState real(pReal), pointer, dimension(:,:) :: & rho, & ! < all dislocations rhoSgl, & - rhoSglMobile, & ! iRhoU + rhoSglMobile, & ! iRhoU rho_sgl_mob_edg_pos, & rho_sgl_mob_edg_neg, & rho_sgl_mob_scr_pos, & @@ -176,14 +176,15 @@ submodule(constitutive) plastic_nonlocal v_scr_pos, & v_scr_neg end type tNonlocalState - + type(tNonlocalState), allocatable, dimension(:) :: & deltaState, & dotState, & - state - + state, & + state0 + type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) - + type(tNonlocalMicrostructure), dimension(:), allocatable :: microstructure contains @@ -210,8 +211,8 @@ module subroutine plastic_nonlocal_init extmsg = '', & structure character(len=pStringLen), dimension(:), allocatable :: outputs - integer :: NofMyPhase - + integer :: NofMyPhase + write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_NONLOCAL_label//' init -+>>>' write(6,'(/,a)') ' Reuber et al., Acta Materialia 71:333–348, 2014' @@ -238,13 +239,14 @@ module subroutine plastic_nonlocal_init associate(prm => param(phase_plasticityInstance(p)), & dot => dotState(phase_plasticityInstance(p)), & stt => state(phase_plasticityInstance(p)), & + st0 => state0(phase_plasticityInstance(p)), & del => deltaState(phase_plasticityInstance(p)), & dst => microstructure(phase_plasticityInstance(p)), & config => config_phase(p)) prm%aTolRho = config%getFloat('atol_rho', defaultVal=0.0_pReal) prm%aTolShear = config%getFloat('atol_shear', defaultVal=0.0_pReal) - + structure = config%getString('lattice_structure') ! This data is read in already in lattice @@ -293,20 +295,20 @@ module subroutine plastic_nonlocal_init prm%colinearSystem(s1) = s2 enddo enddo - + prm%rhoSglEdgePos0 = config%getFloats('rhosgledgepos0', requiredSize=size(prm%Nslip)) prm%rhoSglEdgeNeg0 = config%getFloats('rhosgledgeneg0', requiredSize=size(prm%Nslip)) prm%rhoSglScrewPos0 = config%getFloats('rhosglscrewpos0', requiredSize=size(prm%Nslip)) prm%rhoSglScrewNeg0 = config%getFloats('rhosglscrewneg0', requiredSize=size(prm%Nslip)) prm%rhoDipEdge0 = config%getFloats('rhodipedge0', requiredSize=size(prm%Nslip)) prm%rhoDipScrew0 = config%getFloats('rhodipscrew0', requiredSize=size(prm%Nslip)) - + prm%lambda0 = config%getFloats('lambda0', requiredSize=size(prm%Nslip)) prm%burgers = config%getFloats('burgers', requiredSize=size(prm%Nslip)) - + prm%lambda0 = math_expand(prm%lambda0,prm%Nslip) prm%burgers = math_expand(prm%burgers,prm%Nslip) - + prm%minDipoleHeight_edge = config%getFloats('minimumdipoleheightedge', requiredSize=size(prm%Nslip)) prm%minDipoleHeight_screw = config%getFloats('minimumdipoleheightscrew', requiredSize=size(prm%Nslip)) prm%minDipoleHeight_edge = math_expand(prm%minDipoleHeight_edge,prm%Nslip) @@ -314,7 +316,7 @@ module subroutine plastic_nonlocal_init allocate(prm%minDipoleHeight(prm%totalNslip,2)) prm%minDipoleHeight(:,1) = prm%minDipoleHeight_edge prm%minDipoleHeight(:,2) = prm%minDipoleHeight_screw - + prm%peierlsstress_edge = config%getFloats('peierlsstressedge', requiredSize=size(prm%Nslip)) prm%peierlsstress_screw = config%getFloats('peierlsstressscrew', requiredSize=size(prm%Nslip)) prm%peierlsstress_edge = math_expand(prm%peierlsstress_edge,prm%Nslip) @@ -326,7 +328,7 @@ module subroutine plastic_nonlocal_init prm%significantRho = config%getFloat('significantrho') prm%significantN = config%getFloat('significantn', 0.0_pReal) prm%CFLfactor = config%getFloat('cflfactor',defaultVal=2.0_pReal) - + prm%atomicVolume = config%getFloat('atomicvolume') prm%Dsd0 = config%getFloat('selfdiffusionprefactor') !,'dsd0' prm%selfDiffusionEnergy = config%getFloat('selfdiffusionenergy') !,'qsd' @@ -336,7 +338,7 @@ module subroutine plastic_nonlocal_init prm%solidSolutionEnergy = config%getFloat('solidsolutionenergy') prm%solidSolutionSize = config%getFloat('solidsolutionsize') prm%solidSolutionConcentration = config%getFloat('solidsolutionconcentration') - + prm%p = config%getFloat('p') prm%q = config%getFloat('q') prm%viscosity = config%getFloat('viscosity') @@ -345,62 +347,62 @@ module subroutine plastic_nonlocal_init ! ToDo: discuss logic prm%rhoSglScatter = config%getFloat('rhosglscatter') prm%rhoSglRandom = config%getFloat('rhosglrandom',0.0_pReal) - if (config%keyExists('/rhosglrandom/')) & + if (config%keyExists('/rhosglrandom/')) & prm%rhoSglRandomBinning = config%getFloat('rhosglrandombinning',0.0_pReal) !ToDo: useful default? ! if (rhoSglRandom(instance) < 0.0_pReal) & ! if (rhoSglRandomBinning(instance) <= 0.0_pReal) & - + prm%surfaceTransmissivity = config%getFloat('surfacetransmissivity',defaultVal=1.0_pReal) prm%grainboundaryTransmissivity = config%getFloat('grainboundarytransmissivity',defaultVal=-1.0_pReal) prm%fEdgeMultiplication = config%getFloat('edgemultiplication') prm%shortRangeStressCorrection = config%keyExists('/shortrangestresscorrection/') - + !-------------------------------------------------------------------------------------------------- ! sanity checks if (any(prm%burgers < 0.0_pReal)) extmsg = trim(extmsg)//' burgers' if (any(prm%lambda0 <= 0.0_pReal)) extmsg = trim(extmsg)//' lambda0' - + if (any(prm%rhoSglEdgePos0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoSglEdgePos0' if (any(prm%rhoSglEdgeNeg0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoSglEdgeNeg0' if (any(prm%rhoSglScrewPos0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoSglScrewPos0' if (any(prm%rhoSglScrewNeg0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoSglScrewNeg0' if (any(prm%rhoDipEdge0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoDipEdge0' if (any(prm%rhoDipScrew0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoDipScrew0' - + if (any(prm%peierlsstress < 0.0_pReal)) extmsg = trim(extmsg)//' peierlsstress' if (any(prm%minDipoleHeight < 0.0_pReal)) extmsg = trim(extmsg)//' minDipoleHeight' - + if (prm%viscosity <= 0.0_pReal) extmsg = trim(extmsg)//' viscosity' if (prm%selfDiffusionEnergy <= 0.0_pReal) extmsg = trim(extmsg)//' selfDiffusionEnergy' if (prm%fattack <= 0.0_pReal) extmsg = trim(extmsg)//' fattack' if (prm%doublekinkwidth <= 0.0_pReal) extmsg = trim(extmsg)//' doublekinkwidth' if (prm%Dsd0 < 0.0_pReal) extmsg = trim(extmsg)//' Dsd0' if (prm%atomicVolume <= 0.0_pReal) extmsg = trim(extmsg)//' atomicVolume' ! ToDo: in disloUCLA/dislotwin, the atomic volume is given as a factor - + if (prm%significantN < 0.0_pReal) extmsg = trim(extmsg)//' significantN' if (prm%significantrho < 0.0_pReal) extmsg = trim(extmsg)//' significantrho' if (prm%atolshear <= 0.0_pReal) extmsg = trim(extmsg)//' atolshear' if (prm%atolrho <= 0.0_pReal) extmsg = trim(extmsg)//' atolrho' if (prm%CFLfactor < 0.0_pReal) extmsg = trim(extmsg)//' CFLfactor' - - if (prm%p <= 0.0_pReal .or. prm%p > 1.0_pReal) extmsg = trim(extmsg)//' p' - if (prm%q < 1.0_pReal .or. prm%q > 2.0_pReal) extmsg = trim(extmsg)//' q' - + + if (prm%p <= 0.0_pReal .or. prm%p > 1.0_pReal) extmsg = trim(extmsg)//' p' + if (prm%q < 1.0_pReal .or. prm%q > 2.0_pReal) extmsg = trim(extmsg)//' q' + if (prm%linetensionEffect < 0.0_pReal .or. prm%linetensionEffect > 1.0_pReal) & extmsg = trim(extmsg)//' linetensionEffect' if (prm%edgeJogFactor < 0.0_pReal .or. prm%edgeJogFactor > 1.0_pReal) & extmsg = trim(extmsg)//' edgeJogFactor' - + if (prm%solidSolutionEnergy <= 0.0_pReal) extmsg = trim(extmsg)//' solidSolutionEnergy' if (prm%solidSolutionSize <= 0.0_pReal) extmsg = trim(extmsg)//' solidSolutionSize' if (prm%solidSolutionConcentration <= 0.0_pReal) extmsg = trim(extmsg)//' solidSolutionConcentration' - if (prm%grainboundaryTransmissivity > 1.0_pReal) extmsg = trim(extmsg)//' grainboundaryTransmissivity' + if (prm%grainboundaryTransmissivity > 1.0_pReal) extmsg = trim(extmsg)//' grainboundaryTransmissivity' if (prm%surfaceTransmissivity < 0.0_pReal .or. prm%surfaceTransmissivity > 1.0_pReal) & - extmsg = trim(extmsg)//' surfaceTransmissivity' - + extmsg = trim(extmsg)//' surfaceTransmissivity' + if (prm%fEdgeMultiplication < 0.0_pReal .or. prm%fEdgeMultiplication > 1.0_pReal) & - extmsg = trim(extmsg)//' fEdgeMultiplication' + extmsg = trim(extmsg)//' fEdgeMultiplication' endif slipActive @@ -470,39 +472,40 @@ module subroutine plastic_nonlocal_init 'gamma ' ]) * prm%totalNslip !< "basic" microstructural state variables that are independent from other state variables sizeDependentState = size([ 'rhoForest ']) * prm%totalNslip !< microstructural state variables that depend on other state variables sizeState = sizeDotState + sizeDependentState & - + size([ 'velocityEdgePos ','velocityEdgeNeg ', & + + size([ 'velocityEdgePos ','velocityEdgeNeg ', & 'velocityScrewPos ','velocityScrewNeg ', & 'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%totalNslip !< other dependent state variables that are not updated by microstructure sizeDeltaState = sizeDotState - + call material_allocatePlasticState(p,NofMyPhase,sizeState,sizeDotState,sizeDeltaState) plasticState(p)%nonlocal = .true. plasticState(p)%offsetDeltaState = 0 ! ToDo: state structure does not follow convention - + totalNslip(phase_plasticityInstance(p)) = prm%totalNslip - + + st0%rho => plasticState(p)%state0 (0*prm%totalNslip+1:10*prm%totalNslip,:) stt%rho => plasticState(p)%state (0*prm%totalNslip+1:10*prm%totalNslip,:) dot%rho => plasticState(p)%dotState (0*prm%totalNslip+1:10*prm%totalNslip,:) del%rho => plasticState(p)%deltaState (0*prm%totalNslip+1:10*prm%totalNslip,:) plasticState(p)%aTolState(1:10*prm%totalNslip) = prm%aTolRho - + stt%rhoSgl => plasticState(p)%state (0*prm%totalNslip+1: 8*prm%totalNslip,:) dot%rhoSgl => plasticState(p)%dotState (0*prm%totalNslip+1: 8*prm%totalNslip,:) del%rhoSgl => plasticState(p)%deltaState (0*prm%totalNslip+1: 8*prm%totalNslip,:) - + stt%rhoSglMobile => plasticState(p)%state (0*prm%totalNslip+1: 4*prm%totalNslip,:) dot%rhoSglMobile => plasticState(p)%dotState (0*prm%totalNslip+1: 4*prm%totalNslip,:) del%rhoSglMobile => plasticState(p)%deltaState (0*prm%totalNslip+1: 4*prm%totalNslip,:) - + stt%rho_sgl_mob_edg_pos => plasticState(p)%state (0*prm%totalNslip+1: 1*prm%totalNslip,:) dot%rho_sgl_mob_edg_pos => plasticState(p)%dotState (0*prm%totalNslip+1: 1*prm%totalNslip,:) del%rho_sgl_mob_edg_pos => plasticState(p)%deltaState (0*prm%totalNslip+1: 1*prm%totalNslip,:) - + stt%rho_sgl_mob_edg_neg => plasticState(p)%state (1*prm%totalNslip+1: 2*prm%totalNslip,:) dot%rho_sgl_mob_edg_neg => plasticState(p)%dotState (1*prm%totalNslip+1: 2*prm%totalNslip,:) del%rho_sgl_mob_edg_neg => plasticState(p)%deltaState (1*prm%totalNslip+1: 2*prm%totalNslip,:) - + stt%rho_sgl_mob_scr_pos => plasticState(p)%state (2*prm%totalNslip+1: 3*prm%totalNslip,:) dot%rho_sgl_mob_scr_pos => plasticState(p)%dotState (2*prm%totalNslip+1: 3*prm%totalNslip,:) del%rho_sgl_mob_scr_pos => plasticState(p)%deltaState (2*prm%totalNslip+1: 3*prm%totalNslip,:) @@ -510,46 +513,46 @@ module subroutine plastic_nonlocal_init stt%rho_sgl_mob_scr_neg => plasticState(p)%state (3*prm%totalNslip+1: 4*prm%totalNslip,:) dot%rho_sgl_mob_scr_neg => plasticState(p)%dotState (3*prm%totalNslip+1: 4*prm%totalNslip,:) del%rho_sgl_mob_scr_neg => plasticState(p)%deltaState (3*prm%totalNslip+1: 4*prm%totalNslip,:) - + stt%rhoSglImmobile => plasticState(p)%state (4*prm%totalNslip+1: 8*prm%totalNslip,:) dot%rhoSglImmobile => plasticState(p)%dotState (4*prm%totalNslip+1: 8*prm%totalNslip,:) del%rhoSglImmobile => plasticState(p)%deltaState (4*prm%totalNslip+1: 8*prm%totalNslip,:) - + stt%rho_sgl_imm_edg_pos => plasticState(p)%state (4*prm%totalNslip+1: 5*prm%totalNslip,:) dot%rho_sgl_imm_edg_pos => plasticState(p)%dotState (4*prm%totalNslip+1: 5*prm%totalNslip,:) del%rho_sgl_imm_edg_pos => plasticState(p)%deltaState (4*prm%totalNslip+1: 5*prm%totalNslip,:) - + stt%rho_sgl_imm_edg_neg => plasticState(p)%state (5*prm%totalNslip+1: 6*prm%totalNslip,:) dot%rho_sgl_imm_edg_neg => plasticState(p)%dotState (5*prm%totalNslip+1: 6*prm%totalNslip,:) del%rho_sgl_imm_edg_neg => plasticState(p)%deltaState (5*prm%totalNslip+1: 6*prm%totalNslip,:) - + stt%rho_sgl_imm_scr_pos => plasticState(p)%state (6*prm%totalNslip+1: 7*prm%totalNslip,:) dot%rho_sgl_imm_scr_pos => plasticState(p)%dotState (6*prm%totalNslip+1: 7*prm%totalNslip,:) del%rho_sgl_imm_scr_pos => plasticState(p)%deltaState (6*prm%totalNslip+1: 7*prm%totalNslip,:) - + stt%rho_sgl_imm_scr_neg => plasticState(p)%state (7*prm%totalNslip+1: 8*prm%totalNslip,:) dot%rho_sgl_imm_scr_neg => plasticState(p)%dotState (7*prm%totalNslip+1: 8*prm%totalNslip,:) del%rho_sgl_imm_scr_neg => plasticState(p)%deltaState (7*prm%totalNslip+1: 8*prm%totalNslip,:) - + stt%rhoDip => plasticState(p)%state (8*prm%totalNslip+1:10*prm%totalNslip,:) dot%rhoDip => plasticState(p)%dotState (8*prm%totalNslip+1:10*prm%totalNslip,:) del%rhoDip => plasticState(p)%deltaState (8*prm%totalNslip+1:10*prm%totalNslip,:) - + stt%rho_dip_edg => plasticState(p)%state (8*prm%totalNslip+1: 9*prm%totalNslip,:) dot%rho_dip_edg => plasticState(p)%dotState (8*prm%totalNslip+1: 9*prm%totalNslip,:) del%rho_dip_edg => plasticState(p)%deltaState (8*prm%totalNslip+1: 9*prm%totalNslip,:) - + stt%rho_dip_scr => plasticState(p)%state (9*prm%totalNslip+1:10*prm%totalNslip,:) dot%rho_dip_scr => plasticState(p)%dotState (9*prm%totalNslip+1:10*prm%totalNslip,:) del%rho_dip_scr => plasticState(p)%deltaState (9*prm%totalNslip+1:10*prm%totalNslip,:) - + stt%gamma => plasticState(p)%state (10*prm%totalNslip + 1:11*prm%totalNslip ,1:NofMyPhase) dot%gamma => plasticState(p)%dotState (10*prm%totalNslip + 1:11*prm%totalNslip ,1:NofMyPhase) del%gamma => plasticState(p)%deltaState (10*prm%totalNslip + 1:11*prm%totalNslip ,1:NofMyPhase) plasticState(p)%aTolState(10*prm%totalNslip + 1:11*prm%totalNslip ) = prm%aTolShear plasticState(p)%slipRate => plasticState(p)%dotState (10*prm%totalNslip + 1:11*prm%totalNslip ,1:NofMyPhase) plasticState(p)%accumulatedSlip => plasticState(p)%state(10*prm%totalNslip + 1:11*prm%totalNslip ,1:NofMyPhase) - + stt%rho_forest => plasticState(p)%state (11*prm%totalNslip + 1:12*prm%totalNslip ,1:NofMyPhase) stt%v => plasticState(p)%state (12*prm%totalNslip + 1:16*prm%totalNslip ,1:NofMyPhase) stt%v_edg_pos => plasticState(p)%state (12*prm%totalNslip + 1:13*prm%totalNslip ,1:NofMyPhase) @@ -565,10 +568,10 @@ module subroutine plastic_nonlocal_init plasticState(p)%state0 = plasticState(p)%state enddo - + allocate(compatibility(2,maxval(totalNslip),maxval(totalNslip),nIPneighbors,& discretization_nIP,discretization_nElem), source=0.0_pReal) - + ! BEGIN DEPRECATED---------------------------------------------------------------------------------- allocate(iRhoU(maxval(totalNslip),4,maxNinstances), source=0) allocate(iRhoB(maxval(totalNslip),4,maxNinstances), source=0) @@ -601,7 +604,7 @@ module subroutine plastic_nonlocal_init iRhoD(s,c,phase_plasticityInstance(p)) = l enddo enddo - + l = l + param(phase_plasticityInstance(p))%totalNslip ! shear(rates) l = l + param(phase_plasticityInstance(p))%totalNslip ! rho_forest @@ -619,20 +622,20 @@ module subroutine plastic_nonlocal_init enddo if (iD(param(phase_plasticityInstance(p))%totalNslip,2,phase_plasticityInstance(p)) /= plasticState(p)%sizeState) & call IO_error(0, ext_msg = 'state indices not properly set ('//PLASTICITY_NONLOCAL_label//')') - - + + endif myPhase2 - + enddo initializeInstances ! END DEPRECATED------------------------------------------------------------------------------------ - + contains !-------------------------------------------------------------------------------------------------- !> @brief populates the initial dislocation density !-------------------------------------------------------------------------------------------------- subroutine stateInit(phase,NofMyPhase) - + integer,intent(in) ::& phase, & NofMyPhase @@ -647,7 +650,7 @@ module subroutine plastic_nonlocal_init phasemember real(pReal), dimension(2) :: & noise, & - rnd + rnd real(pReal) :: & meanDensity, & totalVolume, & @@ -655,14 +658,14 @@ module subroutine plastic_nonlocal_init minimumIpVolume real(pReal), dimension(NofMyPhase) :: & volume - - + + instance = phase_plasticityInstance(phase) associate(prm => param(instance), stt => state(instance)) - - ! randomly distribute dislocation segments on random slip system and of random type in the volume + + ! randomly distribute dislocation segments on random slip system and of random type in the volume if (prm%rhoSglRandom > 0.0_pReal) then - + ! get the total volume of the instance do e = 1,discretization_nElem do i = 1,discretization_nIP @@ -672,7 +675,7 @@ module subroutine plastic_nonlocal_init totalVolume = sum(volume) minimumIPVolume = minval(volume) densityBinning = prm%rhoSglRandomBinning / minimumIpVolume ** (2.0_pReal / 3.0_pReal) - + ! subsequently fill random ips with dislocation segments until we reach the desired overall density meanDensity = 0.0_pReal do while(meanDensity < prm%rhoSglRandom) @@ -701,9 +704,9 @@ module subroutine plastic_nonlocal_init enddo enddo endif - + end associate - + end subroutine stateInit end subroutine plastic_nonlocal_init @@ -713,14 +716,14 @@ end subroutine plastic_nonlocal_init !> @brief calculates quantities characterizing the microstructure !-------------------------------------------------------------------------------------------------- module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) - + integer, intent(in) :: & ip, & el real(pReal), dimension(3,3), intent(in) :: & Fe, & Fp - + integer :: & ph, & !< phase of, & !< offset @@ -765,7 +768,7 @@ module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))), & totalNslip(phase_plasticityInstance(material_phaseAt(1,el)))) :: & myInteractionMatrix ! corrected slip interaction matrix - + real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),nIPneighbors) :: & rho_edg_delta_neighbor, & rho_scr_delta_neighbor @@ -774,22 +777,22 @@ module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) neighbor_rhoTotal ! total density at neighboring material point real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),2) :: & m ! direction of dislocation motion - + ph = material_phaseAt(1,el) of = material_phasememberAt(1,ip,el) instance = phase_plasticityInstance(ph) associate(prm => param(instance),dst => microstructure(instance), stt => state(instance)) - + ns = prm%totalNslip - + rho = getRho(instance,of,ip,el) - + stt%rho_forest(:,of) = matmul(prm%forestProjection_Edge, sum(abs(rho(:,edg)),2)) & - + matmul(prm%forestProjection_Screw,sum(abs(rho(:,scr)),2)) - - - ! coefficients are corrected for the line tension effect + + matmul(prm%forestProjection_Screw,sum(abs(rho(:,scr)),2)) + + + ! coefficients are corrected for the line tension effect ! (see Kubin,Devincre,Hoc; 2008; Modeling dislocation storage rates and mean free paths in face-centered cubic crystals) if (lattice_structure(ph) == LATTICE_bcc_ID .or. lattice_structure(ph) == LATTICE_fcc_ID) then ! only fcc and bcc do s = 1,ns @@ -797,7 +800,7 @@ module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) + prm%linetensionEffect & * log(0.35_pReal * prm%burgers(s) * sqrt(max(stt%rho_forest(s,of),prm%significantRho))) & / log(0.35_pReal * prm%burgers(s) * 1e6_pReal)) ** 2.0_pReal - myInteractionMatrix(1:ns,s) = correction * prm%interactionSlipSlip(1:ns,s) + myInteractionMatrix(1:ns,s) = correction * prm%interactionSlipSlip(1:ns,s) enddo else myInteractionMatrix = prm%interactionSlipSlip @@ -819,17 +822,17 @@ module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) if (.not. phase_localPlasticity(ph) .and. prm%shortRangeStressCorrection) then invFe = math_inv33(Fe) invFp = math_inv33(Fp) - + rho_edg_delta = rho(:,mob_edg_pos) - rho(:,mob_edg_neg) rho_scr_delta = rho(:,mob_scr_pos) - rho(:,mob_scr_neg) - + rhoExcess(1,1:ns) = rho_edg_delta rhoExcess(2,1:ns) = rho_scr_delta - + FVsize = IPvolume(ip,el) ** (1.0_pReal/3.0_pReal) - + !* loop through my neighborhood and get the connection vectors (in lattice frame) and the excess densities - + nRealNeighbors = 0.0_pReal neighbor_rhoTotal = 0.0_pReal do n = 1,nIPneighbors @@ -839,16 +842,16 @@ module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) if (neighbor_el > 0 .and. neighbor_ip > 0) then neighbor_instance = phase_plasticityInstance(material_phaseAt(1,neighbor_el)) if (neighbor_instance == instance) then - + nRealNeighbors = nRealNeighbors + 1.0_pReal rho_neighbor = getRho(instance,no,neighbor_ip,neighbor_el) - + rho_edg_delta_neighbor(:,n) = rho_neighbor(:,mob_edg_pos) - rho_neighbor(:,mob_edg_neg) rho_scr_delta_neighbor(:,n) = rho_neighbor(:,mob_scr_pos) - rho_neighbor(:,mob_scr_neg) - + neighbor_rhoTotal(1,:,n) = sum(abs(rho_neighbor(:,edg)),2) neighbor_rhoTotal(2,:,n) = sum(abs(rho_neighbor(:,scr)),2) - + connection_latticeConf(1:3,n) = matmul(invFe, discretization_IPcoords(1:3,neighbor_el+neighbor_ip-1) & - discretization_IPcoords(1:3,el+neighbor_ip-1)) normal_latticeConf = matmul(transpose(invFp), IPareaNormal(1:3,n,ip,el)) @@ -867,18 +870,18 @@ module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) rho_scr_delta_neighbor(:,n) = rho_scr_delta endif enddo - + neighbor_rhoExcess(1,:,:) = rho_edg_delta_neighbor neighbor_rhoExcess(2,:,:) = rho_scr_delta_neighbor - + !* loop through the slip systems and calculate the dislocation gradient by !* 1. interpolation of the excess density in the neighorhood !* 2. interpolation of the dead dislocation density in the central volume m(1:3,1:ns,1) = prm%slip_direction m(1:3,1:ns,2) = -prm%slip_transverse - + do s = 1,ns - + ! gradient from interpolation of neighboring excess density ... do c = 1,2 do dir = 1,3 @@ -891,27 +894,27 @@ module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) enddo invConnections = math_inv33(connections) if (all(dEq0(invConnections))) call IO_error(-1,ext_msg='back stress calculation: inversion error') - + rhoExcessGradient(c) = math_inner(m(1:3,s,c), matmul(invConnections,rhoExcessDifferences)) enddo - + ! ... plus gradient from deads ... rhoExcessGradient(1) = rhoExcessGradient(1) + sum(rho(s,imm_edg)) / FVsize rhoExcessGradient(2) = rhoExcessGradient(2) + sum(rho(s,imm_scr)) / FVsize - + ! ... normalized with the total density ... rhoTotal(1) = (sum(abs(rho(s,edg))) + sum(neighbor_rhoTotal(1,s,:))) / (1.0_pReal + nRealNeighbors) rhoTotal(2) = (sum(abs(rho(s,scr))) + sum(neighbor_rhoTotal(2,s,:))) / (1.0_pReal + nRealNeighbors) - + rhoExcessGradient_over_rho = 0.0_pReal where(rhoTotal > 0.0_pReal) & rhoExcessGradient_over_rho = rhoExcessGradient / rhoTotal - + ! ... gives the local stress correction when multiplied with a factor dst%tau_back(s,of) = - prm%mu * prm%burgers(s) / (2.0_pReal * pi) & * (rhoExcessGradient_over_rho(1) / (1.0_pReal - prm%nu) & + rhoExcessGradient_over_rho(2)) - + enddo endif @@ -932,7 +935,7 @@ end subroutine plastic_nonlocal_dependentState !-------------------------------------------------------------------------------------------------- -!> @brief calculates kinetics +!> @brief calculates kinetics !-------------------------------------------------------------------------------------------------- subroutine plastic_nonlocal_kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, & tauThreshold, c, Temperature, instance, of) @@ -945,17 +948,17 @@ subroutine plastic_nonlocal_kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, & tau, & !< resolved external shear stress (without non Schmid effects) tauNS, & !< resolved external shear stress (including non Schmid effects) tauThreshold !< threshold shear stress - + real(pReal), dimension(param(instance)%totalNslip), intent(out) :: & v, & !< velocity dv_dtau, & !< velocity derivative with respect to resolved shear stress (without non Schmid contributions) dv_dtauNS !< velocity derivative with respect to resolved shear stress (including non Schmid contributions) - + integer :: & ns, & !< short notation for the total number of active slip systems s !< index of my current slip system - real(pReal) :: & - tauRel_P, & + real(pReal) :: & + tauRel_P, & tauRel_S, & tauEff, & !< effective shear stress tPeierls, & !< waiting time in front of a peierls barriers @@ -976,22 +979,22 @@ subroutine plastic_nonlocal_kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, & criticalStress_P, & !< maximum obstacle strength criticalStress_S, & !< maximum obstacle strength mobility !< dislocation mobility - + associate(prm => param(instance)) ns = prm%totalNslip v = 0.0_pReal dv_dtau = 0.0_pReal dv_dtauNS = 0.0_pReal - - + + if (Temperature > 0.0_pReal) then do s = 1,ns if (abs(tau(s)) > tauThreshold(s)) then - + !* Peierls contribution !* Effective stress includes non Schmid constributions !* The derivative only gives absolute values; the correct sign is taken care of in the formula for the derivative of the velocity - + tauEff = max(0.0_pReal, abs(tauNS(s)) - tauThreshold(s)) ! ensure that the effective stress is positive meanfreepath_P = prm%burgers(s) jumpWidth_P = prm%burgers(s) @@ -1006,15 +1009,15 @@ subroutine plastic_nonlocal_kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, & if (tauEff < criticalStress_P) then dtPeierls_dtau = tPeierls * prm%p * prm%q * activationVolume_P / (KB * Temperature) & * (1.0_pReal - tauRel_P**prm%p)**(prm%q-1.0_pReal) & - * tauRel_P**(prm%p-1.0_pReal) + * tauRel_P**(prm%p-1.0_pReal) else dtPeierls_dtau = 0.0_pReal endif - - + + !* Contribution from solid solution strengthening !* The derivative only gives absolute values; the correct sign is taken care of in the formula for the derivative of the velocity - + tauEff = abs(tau(s)) - tauThreshold(s) meanfreepath_S = prm%burgers(s) / sqrt(prm%solidSolutionConcentration) jumpWidth_S = prm%solidSolutionSize * prm%burgers(s) @@ -1030,32 +1033,32 @@ subroutine plastic_nonlocal_kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, & dtSolidSolution_dtau = tSolidSolution * prm%p * prm%q & * activationVolume_S / (KB * Temperature) & * (1.0_pReal - tauRel_S**prm%p)**(prm%q-1.0_pReal) & - * tauRel_S**(prm%p-1.0_pReal) + * tauRel_S**(prm%p-1.0_pReal) else dtSolidSolution_dtau = 0.0_pReal endif - - + + !* viscous glide velocity - + tauEff = abs(tau(s)) - tauThreshold(s) mobility = prm%burgers(s) / prm%viscosity vViscous = mobility * tauEff - - - !* Mean velocity results from waiting time at peierls barriers and solid solution obstacles with respective meanfreepath of - !* free flight at glide velocity in between. + + + !* Mean velocity results from waiting time at peierls barriers and solid solution obstacles with respective meanfreepath of + !* free flight at glide velocity in between. !* adopt sign from resolved stress - + v(s) = sign(1.0_pReal,tau(s)) & / (tPeierls / meanfreepath_P + tSolidSolution / meanfreepath_S + 1.0_pReal / vViscous) dv_dtau(s) = v(s) * v(s) * (dtSolidSolution_dtau / meanfreepath_S & + mobility / (vViscous * vViscous)) - dv_dtauNS(s) = v(s) * v(s) * dtPeierls_dtau / meanfreepath_P + dv_dtauNS(s) = v(s) * v(s) * dtPeierls_dtau / meanfreepath_P endif enddo endif - + #ifdef DEBUGTODO write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> tauThreshold / MPa', tauThreshold * 1e-6_pReal @@ -1089,8 +1092,8 @@ module subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dMp, & Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp !< derivative of Lp with respect to Tstar (9x9 matrix) - - + + integer :: & instance, & !< current instance of this plasticity ns, & !< short notation for the total number of active slip systems @@ -1114,20 +1117,20 @@ module subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dMp, & real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el)))) :: & tau, & !< resolved shear stress including backstress terms gdotTotal !< shear rate - + !*** shortcut for mapping ph = material_phaseAt(1,el) of = material_phasememberAt(1,ip,el) - + instance = phase_plasticityInstance(ph) associate(prm => param(instance),dst=>microstructure(instance),stt=>state(instance)) ns = prm%totalNslip - - !*** shortcut to state variables + + !*** shortcut to state variables rho = getRho(instance,of,ip,el) rhoSgl = rho(:,sgl) - - + + !*** get resolved shear stress !*** for screws possible non-schmid contributions are also taken into account do s = 1,ns @@ -1145,18 +1148,18 @@ module subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dMp, & forall (t = 1:4) & tauNS(1:ns,t) = tauNS(1:ns,t) + dst%tau_back(:,of) tau = tau + dst%tau_back(:,of) - - + + !*** get dislocation velocity and its tangent and store the velocity in the state array - - ! edges + + ! edges call plastic_nonlocal_kinetics(v(1:ns,1), dv_dtau(1:ns,1), dv_dtauNS(1:ns,1), & tau(1:ns), tauNS(1:ns,1), dst%tau_pass(1:ns,of), & 1, Temperature, instance, of) v(1:ns,2) = v(1:ns,1) dv_dtau(1:ns,2) = dv_dtau(1:ns,1) dv_dtauNS(1:ns,2) = dv_dtauNS(1:ns,1) - + !screws if (size(prm%nonSchmidCoeff) == 0) then forall(t = 3:4) @@ -1173,17 +1176,17 @@ module subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dMp, & endif stt%v(:,of) = pack(v,.true.) - + !*** Bauschinger effect forall (s = 1:ns, t = 5:8, rhoSgl(s,t) * v(s,t-4) < 0.0_pReal) & rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t)) - - + + gdotTotal = sum(rhoSgl(1:ns,1:4) * v, 2) * prm%burgers(1:ns) - + Lp = 0.0_pReal dLp_dMp = 0.0_pReal - + do s = 1,ns Lp = Lp + gdotTotal(s) * prm%Schmid(1:3,1:3,s) forall (i=1:3,j=1:3,k=1:3,l=1:3) & @@ -1191,13 +1194,13 @@ module subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dMp, & + prm%Schmid(i,j,s) * prm%Schmid(k,l,s) & * sum(rhoSgl(s,1:4) * dv_dtau(s,1:4)) * prm%burgers(s) & + prm%Schmid(i,j,s) & - * ( prm%nonSchmid_pos(k,l,s) * rhoSgl(s,3) * dv_dtauNS(s,3) & + * ( prm%nonSchmid_pos(k,l,s) * rhoSgl(s,3) * dv_dtauNS(s,3) & - prm%nonSchmid_neg(k,l,s) * rhoSgl(s,4) * dv_dtauNS(s,4)) * prm%burgers(s) enddo - - + + end associate - + end subroutine plastic_nonlocal_LpAndItsTangent @@ -1205,7 +1208,7 @@ end subroutine plastic_nonlocal_LpAndItsTangent !> @brief (instantaneous) incremental change of microstructure !-------------------------------------------------------------------------------------------------- module subroutine plastic_nonlocal_deltaState(Mp,ip,el) - + integer, intent(in) :: & ip, & el @@ -1234,23 +1237,23 @@ module subroutine plastic_nonlocal_deltaState(Mp,ip,el) dUpper, & ! current maximum stable dipole distance for edges and screws dUpperOld, & ! old maximum stable dipole distance for edges and screws deltaDUpper ! change in maximum stable dipole distance for edges and screws - + ph = material_phaseAt(1,el) of = material_phasememberAt(1,ip,el) instance = phase_plasticityInstance(ph) associate(prm => param(instance),dst => microstructure(instance),del => deltaState(instance)) ns = totalNslip(instance) - - !*** shortcut to state variables + + !*** shortcut to state variables forall (s = 1:ns, t = 1:4) & v(s,t) = plasticState(ph)%state(iV(s,t,instance),of) forall (s = 1:ns, c = 1:2) & dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,instance),of) - + rho = getRho(instance,of,ip,el) rhoDip = rho(:,dip) - + !**************************************************************************** !*** dislocation remobilization (bauschinger effect) where(rho(:,imm) * v < 0.0_pReal) @@ -1261,48 +1264,48 @@ module subroutine plastic_nonlocal_deltaState(Mp,ip,el) elsewhere deltaRhoRemobilization(:,mob) = 0.0_pReal deltaRhoRemobilization(:,imm) = 0.0_pReal - endwhere + endwhere deltaRhoRemobilization(:,dip) = 0.0_pReal - + !**************************************************************************** !*** calculate dipole formation and dissociation by stress change - + !*** calculate limits for stable dipole height do s = 1,prm%totalNslip tau(s) = math_mul33xx33(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 - + dUpper(1:ns,1) = prm%mu * prm%burgers/(8.0_pReal * PI * (1.0_pReal - prm%nu) * abs(tau)) dUpper(1:ns,2) = prm%mu * prm%burgers/(4.0_pReal * PI * abs(tau)) - + where(dNeq0(sqrt(sum(abs(rho(:,edg)),2)))) & dUpper(1:ns,1) = min(1.0_pReal/sqrt(sum(abs(rho(:,edg)),2)),dUpper(1:ns,1)) - + where(dNeq0(sqrt(sum(abs(rho(:,scr)),2)))) & dUpper(1:ns,2) = min(1.0_pReal/sqrt(sum(abs(rho(:,scr)),2)),dUpper(1:ns,2)) - + dUpper = max(dUpper,prm%minDipoleHeight) deltaDUpper = dUpper - dUpperOld - - + + !*** dissociation by stress increase deltaRhoDipole2SingleStress = 0.0_pReal forall (c=1:2, s=1:ns, deltaDUpper(s,c) < 0.0_pReal .and. & dNeq0(dUpperOld(s,c) - prm%minDipoleHeight(s,c))) & deltaRhoDipole2SingleStress(s,8+c) = rhoDip(s,c) * deltaDUpper(s,c) & / (dUpperOld(s,c) - prm%minDipoleHeight(s,c)) - + forall (t=1:4) & deltaRhoDipole2SingleStress(1:ns,t) = -0.5_pReal & * deltaRhoDipole2SingleStress(1:ns,(t-1)/2+9) - + forall (s = 1:ns, c = 1:2) & - plasticState(ph)%state(iD(s,c,instance),of) = dUpper(s,c) - + plasticState(ph)%state(iD(s,c,instance),of) = dUpper(s,c) + plasticState(ph)%deltaState(:,of) = 0.0_pReal del%rho(:,of) = reshape(deltaRhoRemobilization + deltaRhoDipole2SingleStress, [10*ns]) - + #ifdef DEBUG if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0 & .and. ((debug_e == el .and. debug_i == ip)& @@ -1322,7 +1325,7 @@ end subroutine plastic_nonlocal_deltaState !--------------------------------------------------------------------------------------------------- module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & timestep,ip,el) - + integer, intent(in) :: & ip, & !< current integration point el !< current element number @@ -1334,9 +1337,9 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: & Fe, & !< elastic deformation gradient Fp !< plastic deformation gradient - + integer :: & - ph, & + ph, & instance, & !< current instance of this plasticity neighbor_instance, & !< instance of my neighbor's plasticity ns, & !< short notation for the total number of active slip systems @@ -1376,7 +1379,7 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el)))) :: & tau, & !< current resolved shear stress vClimb !< climb velocity of edge dipoles - + real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),2) :: & rhoDip, & !< current dipole dislocation densities (screw and edge dipoles) dLower, & !< minimum stable dipole distance for edges and screws @@ -1399,19 +1402,19 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & transmissivity, & !< overall transmissivity of dislocation flux to neighboring material point lineLength, & !< dislocation line length leaving the current interface selfDiffusion !< self diffusion - + logical :: & considerEnteringFlux, & considerLeavingFlux - + p = material_phaseAt(1,el) o = material_phasememberAt(1,ip,el) - + if (timestep <= 0.0_pReal) then plasticState(p)%dotState = 0.0_pReal return endif - + ph = material_phaseAt(1,el) instance = phase_plasticityInstance(ph) associate(prm => param(instance), & @@ -1419,25 +1422,25 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & dot => dotState(instance), & stt => state(instance)) ns = totalNslip(instance) - + tau = 0.0_pReal gdot = 0.0_pReal - + rho = getRho(instance,o,ip,el) rhoSgl = rho(:,sgl) rhoDip = rho(:,dip) - + forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(p)%state(iV (s,t,instance),o) endforall - - + + !**************************************************************************** !*** Calculate shear rate - + forall (t = 1:4) & gdot(1:ns,t) = rhoSgl(1:ns,t) * prm%burgers(1:ns) * v(1:ns,t) - + #ifdef DEBUG if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0 & .and. ((debug_e == el .and. debug_i == ip)& @@ -1446,29 +1449,29 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & write(6,'(a,/,4(12x,12(e12.5,1x),/))') '<< CONST >> gdot / 1/s',gdot endif #endif - - - + + + !**************************************************************************** !*** calculate limits for stable dipole height - + do s = 1,ns ! loop over slip systems tau(s) = math_mul33xx33(Mp, prm%Schmid(1:3,1:3,s)) + dst%tau_back(s,o) if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal enddo - + dLower = prm%minDipoleHeight dUpper(1:ns,1) = prm%mu * prm%burgers/(8.0_pReal * PI * (1.0_pReal - prm%nu) * abs(tau)) dUpper(1:ns,2) = prm%mu * prm%burgers/(4.0_pReal * PI * abs(tau)) - + where(dNeq0(sqrt(sum(abs(rho(:,edg)),2)))) & dUpper(1:ns,1) = min(1.0_pReal/sqrt(sum(abs(rho(:,edg)),2)),dUpper(1:ns,1)) - + where(dNeq0(sqrt(sum(abs(rho(:,scr)),2)))) & dUpper(1:ns,2) = min(1.0_pReal/sqrt(sum(abs(rho(:,scr)),2)),dUpper(1:ns,2)) dUpper = max(dUpper,dLower) - + !**************************************************************************** !*** calculate dislocation multiplication rhoDotMultiplication = 0.0_pReal @@ -1481,25 +1484,25 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & * sqrt(stt%rho_forest(s,o)) / 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:ns,1:4) = spread( & (sum(abs(gdot(1:ns,1:2)),2) * prm%fEdgeMultiplication + sum(abs(gdot(1:ns,3:4)),2)) & * sqrt(stt%rho_forest(:,o)) / prm%lambda0 / prm%burgers(1:ns), 2, 4) endif isBCC - - + + !**************************************************************************** !*** 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(v) * 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 + 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(v), abs(gdot) > 0.0_pReal & @@ -1512,32 +1515,32 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & plasticState(p)%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 p vector in the (s,p,n) triplet !!! - + m(1:3,1:ns,1) = prm%slip_direction m(1:3,1:ns,2) = -prm%slip_direction m(1:3,1:ns,3) = -prm%slip_transverse m(1:3,1:ns,4) = prm%slip_transverse - + my_Fe = Fe(1:3,1:3,1,ip,el) my_F = matmul(my_Fe, 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_Fe = Fe(1:3,1:3,1,neighbor_ip,neighbor_el) @@ -1546,33 +1549,33 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & else ! if no neighbor, take my value as average Favg = my_F endif - - + + !* FLUX FROM MY NEIGHBOR TO ME - !* This is only considered, if I have a neighbor of nonlocal plasticity - !* (also nonlocal constitutive law with local properties) that is at least a little bit + !* This is only considered, if I have a neighbor of nonlocal plasticity + !* (also nonlocal constitutive law with local properties) that is at least a little bit !* compatible. !* If it's not at all compatible, no flux is arriving, because everything is dammed in front of !* my neighbor's interface. - !* The entering flux from my neighbor will be distributed on my slip systems according to the + !* The entering flux from my neighbor will be distributed on my slip systems according to the !* compatibility - + considerEnteringFlux = .false. - neighbor_v = 0.0_pReal ! needed for check of sign change in flux density below + neighbor_v = 0.0_pReal ! needed for check of sign change in flux density below neighbor_rhoSgl = 0.0_pReal if (neighbor_n > 0) then if (phase_plasticity(material_phaseAt(1,neighbor_el)) == PLASTICITY_NONLOCAL_ID & .and. any(compatibility(:,:,:,n,ip,el) > 0.0_pReal)) & considerEnteringFlux = .true. endif - + enteringFlux: if (considerEnteringFlux) then forall (s = 1:ns, t = 1:4) neighbor_v(s,t) = plasticState(np)%state(iV (s,t,neighbor_instance),no) neighbor_rhoSgl(s,t) = max(plasticState(np)%state(iRhoU(s,t,neighbor_instance),no), & 0.0_pReal) endforall - + where (neighbor_rhoSgl * IPvolume(neighbor_ip,neighbor_el) ** 0.667_pReal < prm%significantN & .or. neighbor_rhoSgl < prm%significantRho) & neighbor_rhoSgl = 0.0_pReal @@ -1587,7 +1590,7 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & c = (t + 1) / 2 topp = t + mod(t,2) - mod(t+1,2) if (neighbor_v(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. v(s,t) * neighbor_v(s,t) >= 0.0_pReal ) then ! ... only if no sign change in flux density + .and. v(s,t) * neighbor_v(s,t) >= 0.0_pReal ) then ! ... only if no sign change in flux density lineLength = neighbor_rhoSgl(s,t) * neighbor_v(s,t) & * math_inner(m(1:3,s,t), normal_neighbor2me) * area ! positive line length that wants to enter through this interface where (compatibility(c,1:ns,s,n,ip,el) > 0.0_pReal) & ! positive compatibility... @@ -1602,33 +1605,33 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & enddo enddo endif enteringFlux - - + + !* 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). + !* 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. - + considerLeavingFlux = .true. if (opposite_n > 0) then if (phase_plasticity(material_phaseAt(1,opposite_el)) /= PLASTICITY_NONLOCAL_ID) & considerLeavingFlux = .false. endif - + leavingFlux: if (considerLeavingFlux) then my_rhoSgl = rhoSgl my_v = v - + normal_me2neighbor_defConf = math_det33(Favg) & - * matmul(math_inv33(transpose(Favg)), & + * matmul(math_inv33(transpose(Favg)), & IPareaNormal(1:3,n,ip,el)) ! calculate the normal of the interface in (average) deformed configuration (pointing from me to my 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 + 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 @@ -1642,67 +1645,67 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & * 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, my_v(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 + + lineLength / IPvolume(ip,el) * (1.0_pReal - transmissivity) & + * sign(1.0_pReal, my_v(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 leavingFlux - + enddo neighbors endif - - - + + + !**************************************************************************** !*** calculate dipole formation and annihilation - + !*** formation by glide - + do c = 1,2 rhoDotSingle2DipoleGlide(1:ns,2*c-1) = -2.0_pReal * dUpper(1:ns,c) / prm%burgers(1:ns) & * (rhoSgl(1:ns,2*c-1) * abs(gdot(1:ns,2*c)) & ! negative mobile --> positive mobile + rhoSgl(1:ns,2*c) * abs(gdot(1:ns,2*c-1)) & ! positive mobile --> negative mobile + abs(rhoSgl(1:ns,2*c+4)) * abs(gdot(1:ns,2*c-1))) ! positive mobile --> negative immobile - + rhoDotSingle2DipoleGlide(1:ns,2*c) = -2.0_pReal * dUpper(1:ns,c) / prm%burgers(1:ns) & * (rhoSgl(1:ns,2*c-1) * abs(gdot(1:ns,2*c)) & ! negative mobile --> positive mobile + rhoSgl(1:ns,2*c) * abs(gdot(1:ns,2*c-1)) & ! positive mobile --> negative mobile + abs(rhoSgl(1:ns,2*c+3)) * abs(gdot(1:ns,2*c))) ! negative mobile --> positive immobile - + rhoDotSingle2DipoleGlide(1:ns,2*c+3) = -2.0_pReal * dUpper(1:ns,c) / prm%burgers(1:ns) & * rhoSgl(1:ns,2*c+3) * abs(gdot(1:ns,2*c)) ! negative mobile --> positive immobile - + rhoDotSingle2DipoleGlide(1:ns,2*c+4) = -2.0_pReal * dUpper(1:ns,c) / prm%burgers(1:ns)& * rhoSgl(1:ns,2*c+4) * abs(gdot(1:ns,2*c-1)) ! positive mobile --> negative immobile - + rhoDotSingle2DipoleGlide(1:ns,c+8) = - rhoDotSingle2DipoleGlide(1:ns,2*c-1) & - rhoDotSingle2DipoleGlide(1:ns,2*c) & + abs(rhoDotSingle2DipoleGlide(1:ns,2*c+3)) & + abs(rhoDotSingle2DipoleGlide(1:ns,2*c+4)) enddo - - + + !*** athermal annihilation - + rhoDotAthermalAnnihilation = 0.0_pReal - - forall (c=1:2) & + + forall (c=1:2) & rhoDotAthermalAnnihilation(1:ns,c+8) = -2.0_pReal * dLower(1:ns,c) / prm%burgers(1:ns) & * ( 2.0_pReal * (rhoSgl(1:ns,2*c-1) * abs(gdot(1:ns,2*c)) + rhoSgl(1:ns,2*c) * abs(gdot(1:ns,2*c-1))) & ! was single hitting single + 2.0_pReal * (abs(rhoSgl(1:ns,2*c+3)) * abs(gdot(1:ns,2*c)) + abs(rhoSgl(1:ns,2*c+4)) * abs(gdot(1:ns,2*c-1))) & ! was single hitting immobile single or was immobile single hit by single + rhoDip(1:ns,c) * (abs(gdot(1:ns,2*c-1)) + abs(gdot(1:ns,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,o)) * (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 / ( KB * Temperature ) & @@ -1713,12 +1716,11 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) & - rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have - rhoDot = 0.0_pReal rhoDot = rhoDotFlux & + rhoDotMultiplication & + rhoDotSingle2DipoleGlide & + rhoDotAthermalAnnihilation & - + rhoDotThermalAnnihilation + + rhoDotThermalAnnihilation #ifdef DEBUG if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0 & @@ -1747,7 +1749,7 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & if ( any(rho(:,mob) + rhoDot(1:ns,1:4) * timestep < -prm%aTolRho) & .or. any(rho(:,dip) + rhoDot(1:ns,9:10) * timestep < -prm%aTolRho)) then #ifdef DEBUG - if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0) then + 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 @@ -1765,21 +1767,21 @@ end subroutine plastic_nonlocal_dotState !-------------------------------------------------------------------------------------------------- !> @brief Compatibility update -!> @detail Compatibility is defined as normalized product of signed cosine of the angle between the slip +!> @detail Compatibility is defined as normalized product of signed cosine of the angle between the slip ! plane normals and signed cosine of the angle between the slip directions. Only the largest values ! that sum up to a total of 1 are considered, all others are set to zero. !-------------------------------------------------------------------------------------------------- -module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) - +module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) + integer, intent(in) :: & i, & e type(rotation), dimension(1,discretization_nIP,discretization_nElem), intent(in) :: & orientation ! crystal orientation - + integer :: & Nneighbors, & ! number of neighbors - n, & ! neighbor index + n, & ! neighbor index neighbor_e, & ! element index of my neighbor neighbor_i, & ! integration point index of my neighbor ph, & @@ -1792,13 +1794,13 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) s2 ! slip system index (my neighbor) real(pReal), dimension(2,totalNslip(phase_plasticityInstance(material_phaseAt(1,e))),& totalNslip(phase_plasticityInstance(material_phaseAt(1,e))),& - nIPneighbors) :: & - my_compatibility ! my_compatibility for current element and ip + nIPneighbors) :: & + my_compatibility ! my_compatibility for current element and ip real(pReal) :: & my_compatibilitySum, & thresholdValue, & nThresholdValues - logical, dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,e)))) :: & + logical, dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,e)))) :: & belowThreshold type(rotation) :: mis @@ -1808,32 +1810,32 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) instance = phase_plasticityInstance(ph) ns = totalNslip(instance) associate(prm => param(instance)) - + !*** start out fully compatible my_compatibility = 0.0_pReal - - forall(s1 = 1:ns) my_compatibility(1:2,s1,s1,1:Nneighbors) = 1.0_pReal - + + forall(s1 = 1:ns) my_compatibility(1:2,s1,s1,1:Nneighbors) = 1.0_pReal + !*** Loop thrugh neighbors and check whether there is any compatibility. - + neighbors: do n = 1,Nneighbors neighbor_e = IPneighborhood(1,n,i,e) neighbor_i = IPneighborhood(2,n,i,e) - - + + !* FREE SURFACE !* Set surface transmissivity to the value specified in the material.config - + if (neighbor_e <= 0 .or. neighbor_i <= 0) then forall(s1 = 1:ns) my_compatibility(1:2,s1,s1,n) = sqrt(prm%surfaceTransmissivity) cycle endif - - + + !* PHASE BOUNDARY - !* If we encounter a different nonlocal phase at the neighbor, + !* If we encounter a different nonlocal phase at the neighbor, !* we consider this to be a real "physical" phase boundary, so completely incompatible. - !* If one of the two phases has a local plasticity law, + !* If one of the two phases has a local plasticity law, !* we do not consider this to be a phase boundary, so completely compatible. neighbor_phase = material_phaseAt(1,neighbor_e) if (neighbor_phase /= ph) then @@ -1841,8 +1843,8 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) forall(s1 = 1:ns) my_compatibility(1:2,s1,s1,n) = 0.0_pReal cycle endif - - + + !* GRAIN BOUNDARY ! !* fixed transmissivity for adjacent ips with different texture (only if explicitly given in material.config) if (prm%grainboundaryTransmissivity >= 0.0_pReal) then @@ -1854,16 +1856,16 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) endif cycle endif - - + + !* GRAIN BOUNDARY ? !* Compatibility defined by relative orientation of slip systems: !* The my_compatibility value is defined as the product of the slip normal projection and the slip direction projection. - !* Its sign is always positive for screws, for edges it has the same sign as the slip normal projection. - !* Since the sum for each slip system can easily exceed one (which would result in a transmissivity larger than one), + !* Its sign is always positive for screws, for edges it has the same sign as the slip normal projection. + !* Since the sum for each slip system can easily exceed one (which would result in a transmissivity larger than one), !* only values above or equal to a certain threshold value are considered. This threshold value is chosen, such that - !* the number of compatible slip systems is minimized with the sum of the original compatibility values exceeding one. - !* Finally the smallest compatibility value is decreased until the sum is exactly equal to one. + !* the number of compatible slip systems is minimized with the sum of the original compatibility values exceeding one. + !* Finally the smallest compatibility value is decreased until the sum is exactly equal to one. !* All values below the threshold are set to zero. else mis = orientation(1,i,e)%misorientation(orientation(1,neighbor_i,neighbor_e)) @@ -1878,7 +1880,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) * abs(math_inner(prm%slip_direction(1:3,s1), & mis%rotate(prm%slip_direction(1:3,s2)))) enddo neighborSlipSystems - + my_compatibilitySum = 0.0_pReal belowThreshold = .true. do while (my_compatibilitySum < 1.0_pReal .and. any(belowThreshold(1:ns))) @@ -1896,33 +1898,33 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) where (belowThreshold(1:ns)) my_compatibility(2,1:ns,s1,n) = 0.0_pReal enddo mySlipSystems endif - + enddo neighbors - + compatibility(1:2,1:ns,1:ns,1:Nneighbors,i,e) = my_compatibility - + end associate - + end subroutine plastic_nonlocal_updateCompatibility - + !-------------------------------------------------------------------------------------------------- !> @brief returns copy of current dislocation densities from state !> @details raw values is rectified !-------------------------------------------------------------------------------------------------- function getRho(instance,of,ip,el) - + integer, intent(in) :: instance, of,ip,el real(pReal), dimension(param(instance)%totalNslip,10) :: getRho - + associate(prm => param(instance)) getRho = reshape(state(instance)%rho(:,of),[prm%totalNslip,10]) - + ! ensure positive densities (not for imm, they have a sign) getRho(:,mob) = max(getRho(:,mob),0.0_pReal) getRho(:,dip) = max(getRho(:,dip),0.0_pReal) - + where(abs(getRho) < max(prm%significantN/IPvolume(ip,el)**(2.0_pReal/3.0_pReal),prm%significantRho)) & getRho = 0.0_pReal @@ -1931,6 +1933,31 @@ function getRho(instance,of,ip,el) end function getRho +!-------------------------------------------------------------------------------------------------- +!> @brief returns copy of current dislocation densities from state +!> @details raw values is rectified +!-------------------------------------------------------------------------------------------------- +function getRho0(instance,of,ip,el) + + integer, intent(in) :: instance, of,ip,el + real(pReal), dimension(param(instance)%totalNslip,10) :: getRho0 + + associate(prm => param(instance)) + + getRho0 = reshape(state0(instance)%rho(:,of),[prm%totalNslip,10]) + + ! ensure positive densities (not for imm, they have a sign) + getRho0(:,mob) = max(getRho0(:,mob),0.0_pReal) + getRho0(:,dip) = max(getRho0(:,dip),0.0_pReal) + + where(abs(getRho0) < max(prm%significantN/IPvolume(ip,el)**(2.0_pReal/3.0_pReal),prm%significantRho)) & + getRho0 = 0.0_pReal + + end associate + +end function getRho0 + + !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- From 4f4c6c59496c80d162bf4e0b586a6120ef05fd64 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 7 Feb 2020 11:53:50 +0100 Subject: [PATCH 2/7] using converged dislocation velocity (from last subinc) --- src/constitutive_plastic_nonlocal.f90 | 33 +++++++++++++++------------ 1 file changed, 18 insertions(+), 15 deletions(-) diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index 358fd964d..b7bf340fb 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -227,6 +227,7 @@ module subroutine plastic_nonlocal_init allocate(param(maxNinstances)) allocate(state(maxNinstances)) + allocate(state0(maxNinstances)) allocate(dotState(maxNinstances)) allocate(deltaState(maxNinstances)) allocate(microstructure(maxNinstances)) @@ -1373,8 +1374,8 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & my_rhoSgl !< single dislocation densities of central ip (positive/negative screw and edge without dipoles) real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),4) :: & v, & !< current dislocation glide velocity - my_v, & !< dislocation glide velocity of central ip - neighbor_v, & !< dislocation glide velocity of enighboring ip + v0, & + neighbor_v0, & !< dislocation glide velocity of enighboring ip gdot !< shear rates real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el)))) :: & tau, & !< current resolved shear stress @@ -1491,6 +1492,9 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & * sqrt(stt%rho_forest(:,o)) / prm%lambda0 / prm%burgers(1:ns), 2, 4) endif isBCC + forall (s = 1:ns, t = 1:4) + v0(s,t) = plasticState(p)%state0(iV(s,t,instance),o) + endforall !**************************************************************************** !*** calculate dislocation fluxes (only for nonlocal plasticity) @@ -1499,14 +1503,14 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & !*** check CFL (Courant-Friedrichs-Lewy) condition for flux if (any( abs(gdot) > 0.0_pReal & ! any active slip system ... - .and. prm%CFLfactor * abs(v) * timestep & + .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(v), abs(gdot) > 0.0_pReal & - .and. prm%CFLfactor * abs(v) * timestep & + 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 !!!' @@ -1561,7 +1565,7 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & !* compatibility considerEnteringFlux = .false. - neighbor_v = 0.0_pReal ! needed for check of sign change in flux density below + neighbor_v0 = 0.0_pReal ! needed for check of sign change in flux density below neighbor_rhoSgl = 0.0_pReal if (neighbor_n > 0) then if (phase_plasticity(material_phaseAt(1,neighbor_el)) == PLASTICITY_NONLOCAL_ID & @@ -1571,7 +1575,7 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & enteringFlux: if (considerEnteringFlux) then forall (s = 1:ns, t = 1:4) - neighbor_v(s,t) = plasticState(np)%state(iV (s,t,neighbor_instance),no) + neighbor_v0(s,t) = plasticState(np)%state0(iV (s,t,neighbor_instance),no) neighbor_rhoSgl(s,t) = max(plasticState(np)%state(iRhoU(s,t,neighbor_instance),no), & 0.0_pReal) endforall @@ -1589,9 +1593,9 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & do t = 1,4 c = (t + 1) / 2 topp = t + mod(t,2) - mod(t+1,2) - if (neighbor_v(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. v(s,t) * neighbor_v(s,t) >= 0.0_pReal ) then ! ... only if no sign change in flux density - lineLength = neighbor_rhoSgl(s,t) * neighbor_v(s,t) & + 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_rhoSgl(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,1:ns,s,n,ip,el) > 0.0_pReal) & ! positive compatibility... rhoDotFlux(1:ns,t) = rhoDotFlux(1:ns,t) & @@ -1623,7 +1627,6 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & leavingFlux: if (considerLeavingFlux) then my_rhoSgl = rhoSgl - my_v = v normal_me2neighbor_defConf = math_det33(Favg) & * matmul(math_inv33(transpose(Favg)), & @@ -1635,18 +1638,18 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & do s = 1,ns do t = 1,4 c = (t + 1) / 2 - if (my_v(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 (my_v(s,t) * neighbor_v(s,t) >= 0.0_pReal) then ! no sign change in flux density + 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,1:ns,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_rhoSgl(s,t) * my_v(s,t) & + lineLength = my_rhoSgl(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, my_v(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 + * 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 From f854dc27e941b2a180650abf3325d9c3deaacf03 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 7 Feb 2020 12:23:22 +0100 Subject: [PATCH 3/7] explicit dotState for nonlocal all flux related quantities are calculated based on the converged quantities --- src/constitutive.f90 | 194 ++++++++++----------- src/constitutive_plastic_nonlocal.f90 | 25 ++- src/crystallite.f90 | 234 +++++++++++++------------- 3 files changed, 226 insertions(+), 227 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 681691c9e..f89c08be8 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -25,44 +25,44 @@ module constitutive use kinematics_cleavage_opening use kinematics_slipplane_opening use kinematics_thermal_expansion - + implicit none private - + integer, public, protected :: & constitutive_plasticity_maxSizeDotState, & constitutive_source_maxSizeDotState - + interface module subroutine plastic_none_init end subroutine plastic_none_init - + module subroutine plastic_isotropic_init end subroutine plastic_isotropic_init - + module subroutine plastic_phenopowerlaw_init end subroutine plastic_phenopowerlaw_init - + module subroutine plastic_kinehardening_init end subroutine plastic_kinehardening_init - + module subroutine plastic_dislotwin_init end subroutine plastic_dislotwin_init module subroutine plastic_disloUCLA_init end subroutine plastic_disloUCLA_init - + module subroutine plastic_nonlocal_init end subroutine plastic_nonlocal_init - - + + module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp !< derivative of Lp with respect to the Mandel stress - + real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & @@ -75,20 +75,20 @@ module constitutive Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp !< derivative of Lp with respect to the Mandel stress - + real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & instance, & of end subroutine plastic_phenopowerlaw_LpAndItsTangent - + pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp !< derivative of Lp with respect to the Mandel stress - + real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & @@ -101,7 +101,7 @@ module constitutive Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp !< derivative of Lp with respect to the Mandel stress - + real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress real(pReal), intent(in) :: & @@ -110,13 +110,13 @@ module constitutive instance, & of end subroutine plastic_dislotwin_LpAndItsTangent - + pure module subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp !< derivative of Lp with respect to the Mandel stress - + real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress real(pReal), intent(in) :: & @@ -125,14 +125,14 @@ module constitutive instance, & of end subroutine plastic_disloUCLA_LpAndItsTangent - + module subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dMp, & Mp, Temperature, volume, ip, el) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp !< derivative of Lp with respect to the Mandel stress - + real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress real(pReal), intent(in) :: & @@ -149,15 +149,15 @@ module constitutive Li !< inleastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLi_dMi !< derivative of Li with respect to Mandel stress - + real(pReal), dimension(3,3), intent(in) :: & - Mi !< Mandel stress + Mi !< Mandel stress integer, intent(in) :: & instance, & of end subroutine plastic_isotropic_LiAndItsTangent - - + + module subroutine plastic_isotropic_dotState(Mp,instance,of) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -165,7 +165,7 @@ module constitutive instance, & of end subroutine plastic_isotropic_dotState - + module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -173,7 +173,7 @@ module constitutive instance, & of end subroutine plastic_phenopowerlaw_dotState - + module subroutine plastic_kinehardening_dotState(Mp,instance,of) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -181,7 +181,7 @@ module constitutive instance, & of end subroutine plastic_kinehardening_dotState - + module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -201,8 +201,8 @@ module constitutive instance, & of end subroutine plastic_disloUCLA_dotState - - module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & + + module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature, & timestep,ip,el) integer, intent(in) :: & ip, & !< current integration point @@ -213,11 +213,11 @@ module constitutive real(pReal), dimension(3,3), intent(in) ::& Mp !< MandelStress real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: & - Fe, & !< elastic deformation gradient + F, & !< deformation gradient Fp !< plastic deformation gradient end subroutine plastic_nonlocal_dotState - + module subroutine plastic_dislotwin_dependentState(T,instance,of) integer, intent(in) :: & instance, & @@ -225,13 +225,13 @@ module constitutive real(pReal), intent(in) :: & T end subroutine plastic_dislotwin_dependentState - + module subroutine plastic_disloUCLA_dependentState(instance,of) integer, intent(in) :: & instance, & of end subroutine plastic_disloUCLA_dependentState - + module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) integer, intent(in) :: & ip, & @@ -249,7 +249,7 @@ module constitutive instance, & of end subroutine plastic_kinehardening_deltaState - + module subroutine plastic_nonlocal_deltaState(Mp,ip,el) integer, intent(in) :: & ip, & @@ -267,48 +267,48 @@ module constitutive ip, & !< integration point el !< element end function plastic_dislotwin_homogenizedC - - module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) + + module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) integer, intent(in) :: & i, & e type(rotation), dimension(1,discretization_nIP,discretization_nElem), intent(in) :: & orientation !< crystal orientation end subroutine plastic_nonlocal_updateCompatibility - - + + module subroutine plastic_isotropic_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group end subroutine plastic_isotropic_results - + module subroutine plastic_phenopowerlaw_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group end subroutine plastic_phenopowerlaw_results - + module subroutine plastic_kinehardening_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group end subroutine plastic_kinehardening_results - + module subroutine plastic_dislotwin_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group end subroutine plastic_dislotwin_results - + module subroutine plastic_disloUCLA_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group end subroutine plastic_disloUCLA_results - + module subroutine plastic_nonlocal_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group end subroutine plastic_nonlocal_results - + end interface - + public :: & plastic_nonlocal_updateCompatibility, & constitutive_init, & @@ -321,7 +321,7 @@ module constitutive constitutive_collectDotState, & constitutive_collectDeltaState, & constitutive_results - + contains @@ -355,18 +355,18 @@ subroutine constitutive_init if (any(phase_source == SOURCE_damage_isoDuctile_ID)) call source_damage_isoDuctile_init if (any(phase_source == SOURCE_damage_anisoBrittle_ID)) call source_damage_anisoBrittle_init if (any(phase_source == SOURCE_damage_anisoDuctile_ID)) call source_damage_anisoDuctile_init - + !-------------------------------------------------------------------------------------------------- ! initialize kinematic mechanisms if (any(phase_kinematics == KINEMATICS_cleavage_opening_ID)) call kinematics_cleavage_opening_init if (any(phase_kinematics == KINEMATICS_slipplane_opening_ID)) call kinematics_slipplane_opening_init if (any(phase_kinematics == KINEMATICS_thermal_expansion_ID)) call kinematics_thermal_expansion_init - + write(6,'(/,a)') ' <<<+- constitutive init -+>>>'; flush(6) - + constitutive_plasticity_maxSizeDotState = 0 constitutive_source_maxSizeDotState = 0 - + PhaseLoop2:do ph = 1,material_Nphase !-------------------------------------------------------------------------------------------------- ! partition and inititalize state @@ -398,7 +398,7 @@ function constitutive_homogenizedC(ipc,ip,el) ipc, & !< component-ID of integration point ip, & !< integration point el !< element - + plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) case (PLASTICITY_DISLOTWIN_ID) plasticityType constitutive_homogenizedC = plastic_dislotwin_homogenizedC(ipc,ip,el) @@ -425,10 +425,10 @@ subroutine constitutive_microstructure(Fe, Fp, ipc, ip, el) ho, & !< homogenization tme, & !< thermal member position instance, of - + ho = material_homogenizationAt(el) tme = thermalMapping(ho)%p(ip,el) - + plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) case (PLASTICITY_DISLOTWIN_ID) plasticityType of = material_phasememberAt(ipc,ip,el) @@ -451,7 +451,7 @@ end subroutine constitutive_microstructure ! Mp in, dLp_dMp out !-------------------------------------------------------------------------------------------------- subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & - S, Fi, ipc, ip, el) + S, Fi, ipc, ip, el) integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point @@ -473,49 +473,49 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & tme !< thermal member position integer :: & i, j, instance, of - + ho = material_homogenizationAt(el) tme = thermalMapping(ho)%p(ip,el) - + Mp = matmul(matmul(transpose(Fi),Fi),S) - + plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) - + case (PLASTICITY_NONE_ID) plasticityType Lp = 0.0_pReal dLp_dMp = 0.0_pReal - + case (PLASTICITY_ISOTROPIC_ID) plasticityType of = material_phasememberAt(ipc,ip,el) instance = phase_plasticityInstance(material_phaseAt(ipc,el)) call plastic_isotropic_LpAndItsTangent (Lp,dLp_dMp,Mp,instance,of) - + case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType of = material_phasememberAt(ipc,ip,el) instance = phase_plasticityInstance(material_phaseAt(ipc,el)) call plastic_phenopowerlaw_LpAndItsTangent (Lp,dLp_dMp,Mp,instance,of) - + case (PLASTICITY_KINEHARDENING_ID) plasticityType of = material_phasememberAt(ipc,ip,el) instance = phase_plasticityInstance(material_phaseAt(ipc,el)) call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMp, Mp,instance,of) - + case (PLASTICITY_NONLOCAL_ID) plasticityType call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMp,Mp, & temperature(ho)%p(tme),geometry_plastic_nonlocal_IPvolume0(ip,el),ip,el) - + case (PLASTICITY_DISLOTWIN_ID) plasticityType of = material_phasememberAt(ipc,ip,el) instance = phase_plasticityInstance(material_phaseAt(ipc,el)) call plastic_dislotwin_LpAndItsTangent (Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,of) - + case (PLASTICITY_DISLOUCLA_ID) plasticityType of = material_phasememberAt(ipc,ip,el) instance = phase_plasticityInstance(material_phaseAt(ipc,el)) call plastic_disloucla_LpAndItsTangent (Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,of) - + end select plasticityType - + do i=1,3; do j=1,3 dLp_dFi(i,j,1:3,1:3) = matmul(matmul(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + & matmul(matmul(Fi,dLp_dMp(i,j,1:3,1:3)),S) @@ -545,7 +545,7 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & real(pReal), intent(out), dimension(3,3,3,3) :: & dLi_dS, & !< derivative of Li with respect to S dLi_dFi - + real(pReal), dimension(3,3) :: & my_Li, & !< intermediate velocity gradient FiInv, & @@ -557,11 +557,11 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & integer :: & k, i, j, & instance, of - + Li = 0.0_pReal dLi_dS = 0.0_pReal dLi_dFi = 0.0_pReal - + plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) case (PLASTICITY_isotropic_ID) plasticityType of = material_phasememberAt(ipc,ip,el) @@ -571,10 +571,10 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & my_Li = 0.0_pReal my_dLi_dS = 0.0_pReal end select plasticityType - + Li = Li + my_Li dLi_dS = dLi_dS + my_dLi_dS - + KinematicsLoop: do k = 1, phase_Nkinematics(material_phaseAt(ipc,el)) kinematicsType: select case (phase_kinematics(k,material_phaseAt(ipc,el))) case (KINEMATICS_cleavage_opening_ID) kinematicsType @@ -590,12 +590,12 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & Li = Li + my_Li dLi_dS = dLi_dS + my_dLi_dS enddo KinematicsLoop - + FiInv = math_inv33(Fi) detFi = math_det33(Fi) Li = matmul(matmul(Fi,Li),FiInv)*detFi !< push forward to intermediate configuration temp_33 = matmul(FiInv,Li) - + do i = 1,3; do j = 1,3 dLi_dS(1:3,1:3,i,j) = matmul(matmul(Fi,dLi_dS(1:3,1:3,i,j)),FiInv)*detFi dLi_dFi(1:3,1:3,i,j) = dLi_dFi(1:3,1:3,i,j) + Li*FiInv(j,i) @@ -621,10 +621,10 @@ pure function constitutive_initialFi(ipc, ip, el) integer :: & phase, & homog, offset - + constitutive_initialFi = math_I3 phase = material_phaseAt(ipc,el) - + KinematicsLoop: do k = 1, phase_Nkinematics(phase) !< Warning: small initial strain assumption kinematicsType: select case (phase_kinematics(k,phase)) case (KINEMATICS_thermal_expansion_ID) kinematicsType @@ -640,7 +640,7 @@ end function constitutive_initialFi !-------------------------------------------------------------------------------------------------- !> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to -!> the elastic/intermediate deformation gradients depending on the selected elastic law +!> the elastic/intermediate deformation gradients depending on the selected elastic law !! (so far no case switch because only Hooke is implemented) !-------------------------------------------------------------------------------------------------- subroutine constitutive_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip, el) @@ -690,17 +690,17 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & d !< counter in degradation loop integer :: & i, j - + ho = material_homogenizationAt(el) C = math_66toSym3333(constitutive_homogenizedC(ipc,ip,el)) - + DegradationLoop: do d = 1, phase_NstiffnessDegradations(material_phaseAt(ipc,el)) degradationType: select case(phase_stiffnessDegradation(d,material_phaseAt(ipc,el))) case (STIFFNESS_DEGRADATION_damage_ID) degradationType C = C * damage(ho)%p(damageMapping(ho)%p(ip,el))**2 end select degradationType enddo DegradationLoop - + E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration S = math_mul3333xx33(C,matmul(matmul(transpose(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration @@ -715,7 +715,7 @@ end subroutine constitutive_hooke_SandItsTangents !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip, el) +subroutine constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el) integer, intent(in) :: & ipc, & !< component-ID of integration point @@ -724,7 +724,7 @@ subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip, real(pReal), intent(in) :: & subdt !< timestep real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & - FeArray, & !< elastic deformation gradient + FArray, & !< elastic deformation gradient FpArray !< plastic deformation gradient real(pReal), intent(in), dimension(3,3) :: & Fi !< intermediate deformation gradient @@ -771,7 +771,7 @@ subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip, call plastic_disloucla_dotState (Mp,temperature(ho)%p(tme),instance,of) case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_dotState (Mp,FeArray,FpArray,temperature(ho)%p(tme), & + call plastic_nonlocal_dotState (Mp,FArray,FpArray,temperature(ho)%p(tme), & subdt,ip,el) end select plasticityType @@ -858,32 +858,32 @@ subroutine constitutive_results do p=1,size(config_name_phase) group = trim('current/constituent')//'/'//trim(config_name_phase(p)) call results_closeGroup(results_addGroup(group)) - + group = trim(group)//'/plastic' - - call results_closeGroup(results_addGroup(group)) + + call results_closeGroup(results_addGroup(group)) select case(phase_plasticity(p)) - + case(PLASTICITY_ISOTROPIC_ID) - call plastic_isotropic_results(phase_plasticityInstance(p),group) - + call plastic_isotropic_results(phase_plasticityInstance(p),group) + case(PLASTICITY_PHENOPOWERLAW_ID) - call plastic_phenopowerlaw_results(phase_plasticityInstance(p),group) - + call plastic_phenopowerlaw_results(phase_plasticityInstance(p),group) + case(PLASTICITY_KINEHARDENING_ID) - call plastic_kinehardening_results(phase_plasticityInstance(p),group) + call plastic_kinehardening_results(phase_plasticityInstance(p),group) case(PLASTICITY_DISLOTWIN_ID) - call plastic_dislotwin_results(phase_plasticityInstance(p),group) - + call plastic_dislotwin_results(phase_plasticityInstance(p),group) + case(PLASTICITY_DISLOUCLA_ID) - call plastic_disloUCLA_results(phase_plasticityInstance(p),group) - + call plastic_disloUCLA_results(phase_plasticityInstance(p),group) + case(PLASTICITY_NONLOCAL_ID) - call plastic_nonlocal_results(phase_plasticityInstance(p),group) + call plastic_nonlocal_results(phase_plasticityInstance(p),group) end select - - enddo + + enddo end subroutine constitutive_results diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index b7bf340fb..0a51e78d9 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -1324,7 +1324,7 @@ end subroutine plastic_nonlocal_deltaState !--------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !--------------------------------------------------------------------------------------------------- -module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & +module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature, & timestep,ip,el) integer, intent(in) :: & @@ -1336,7 +1336,7 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & real(pReal), dimension(3,3), intent(in) ::& Mp !< MandelStress real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: & - Fe, & !< elastic deformation gradient + F, & !< elastic deformation gradient Fp !< plastic deformation gradient integer :: & @@ -1370,7 +1370,7 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & rhoDotThermalAnnihilation !< density evolution by thermal annihilation real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),8) :: & rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles) - neighbor_rhoSgl, & !< current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles) + neighbor_rhoSgl0, & !< current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles) my_rhoSgl !< single dislocation densities of central ip (positive/negative screw and edge without dipoles) real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),4) :: & v, & !< current dislocation glide velocity @@ -1529,8 +1529,8 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & m(1:3,1:ns,3) = -prm%slip_transverse m(1:3,1:ns,4) = prm%slip_transverse - my_Fe = Fe(1:3,1:3,1,ip,el) - my_F = matmul(my_Fe, Fp(1:3,1:3,1,ip,el)) + 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 @@ -1547,8 +1547,8 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient neighbor_instance = phase_plasticityInstance(material_phaseAt(1,neighbor_el)) - neighbor_Fe = Fe(1:3,1:3,1,neighbor_ip,neighbor_el) - neighbor_F = matmul(neighbor_Fe, Fp(1:3,1:3,1,neighbor_ip,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 @@ -1566,7 +1566,6 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & considerEnteringFlux = .false. neighbor_v0 = 0.0_pReal ! needed for check of sign change in flux density below - neighbor_rhoSgl = 0.0_pReal 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)) & @@ -1576,13 +1575,13 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & enteringFlux: if (considerEnteringFlux) then forall (s = 1:ns, t = 1:4) neighbor_v0(s,t) = plasticState(np)%state0(iV (s,t,neighbor_instance),no) - neighbor_rhoSgl(s,t) = max(plasticState(np)%state(iRhoU(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_rhoSgl * IPvolume(neighbor_ip,neighbor_el) ** 0.667_pReal < prm%significantN & - .or. neighbor_rhoSgl < prm%significantRho) & - neighbor_rhoSgl = 0.0_pReal + 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)) ! calculate the normal of the interface in (average) deformed configuration (now pointing from my neighbor to me!!!) normal_neighbor2me = matmul(transpose(neighbor_Fe), normal_neighbor2me_defConf) & @@ -1595,7 +1594,7 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & 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_rhoSgl(s,t) * neighbor_v0(s,t) & + 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,1:ns,s,n,ip,el) > 0.0_pReal) & ! positive compatibility... rhoDotFlux(1:ns,t) = rhoDotFlux(1:ns,t) & diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 8356d0a65..53164c8ee 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -22,10 +22,10 @@ module crystallite use discretization use lattice use results - - implicit none + + implicit none private - + real(pReal), dimension(:,:,:), allocatable, public :: & crystallite_dt !< requested time increment of each grain real(pReal), dimension(:,:,:), allocatable :: & @@ -33,7 +33,7 @@ module crystallite crystallite_subFrac, & !< already calculated fraction of increment crystallite_subStep !< size of next integration step type(rotation), dimension(:,:,:), allocatable :: & - crystallite_orientation !< current orientation + crystallite_orientation !< current orientation real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: & crystallite_Fe, & !< current "elastic" def grad (end of converged time step) crystallite_P !< 1st Piola-Kirchhoff stress per grain @@ -74,33 +74,33 @@ module crystallite crystallite_converged, & !< convergence flag crystallite_todo, & !< flag to indicate need for further computation crystallite_localPlasticity !< indicates this grain to have purely local constitutive law - + type :: tOutput !< new requested output (per phase) character(len=pStringLen), allocatable, dimension(:) :: & label end type tOutput type(tOutput), allocatable, dimension(:) :: output_constituent - + type :: tNumerics integer :: & iJacoLpresiduum, & !< frequency of Jacobian update of residuum in Lp nState, & !< state loop limit - nStress !< stress loop limit + nStress !< stress loop limit real(pReal) :: & subStepMinCryst, & !< minimum (relative) size of sub-step allowed during cutback subStepSizeCryst, & !< size of first substep when cutback subStepSizeLp, & !< size of first substep when cutback in Lp calculation subStepSizeLi, & !< size of first substep when cutback in Li calculation stepIncreaseCryst, & !< increase of next substep size when previous substep converged - rTol_crystalliteState, & !< relative tolerance in state loop + rTol_crystalliteState, & !< relative tolerance in state loop rTol_crystalliteStress, & !< relative tolerance in stress loop aTol_crystalliteStress !< absolute tolerance in stress loop end type tNumerics - + type(tNumerics) :: num ! numerics parameters. Better name? - + procedure(), pointer :: integrateState - + public :: & crystallite_init, & crystallite_stress, & @@ -116,7 +116,7 @@ contains !> @brief allocates and initialize per grain variables !-------------------------------------------------------------------------------------------------- subroutine crystallite_init - + logical, dimension(discretization_nIP,discretization_nElem) :: devNull integer :: & c, & !< counter in integration point component loop @@ -126,13 +126,13 @@ subroutine crystallite_init iMax, & !< maximum number of integration points eMax, & !< maximum number of elements myNcomponents !< number of components at current IP - + write(6,'(/,a)') ' <<<+- crystallite init -+>>>' - + cMax = homogenization_maxNgrains iMax = discretization_nIP eMax = discretization_nElem - + allocate(crystallite_S0(3,3,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_partionedS0(3,3,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_S(3,3,cMax,iMax,eMax), source=0.0_pReal) @@ -172,23 +172,23 @@ subroutine crystallite_init allocate(crystallite_requested(cMax,iMax,eMax), source=.false.) allocate(crystallite_todo(cMax,iMax,eMax), source=.false.) allocate(crystallite_converged(cMax,iMax,eMax), source=.true.) - + num%subStepMinCryst = config_numerics%getFloat('substepmincryst', defaultVal=1.0e-3_pReal) num%subStepSizeCryst = config_numerics%getFloat('substepsizecryst', defaultVal=0.25_pReal) num%stepIncreaseCryst = config_numerics%getFloat('stepincreasecryst', defaultVal=1.5_pReal) - + num%subStepSizeLp = config_numerics%getFloat('substepsizelp', defaultVal=0.5_pReal) num%subStepSizeLi = config_numerics%getFloat('substepsizeli', defaultVal=0.5_pReal) num%rTol_crystalliteState = config_numerics%getFloat('rtol_crystallitestate', defaultVal=1.0e-6_pReal) num%rTol_crystalliteStress = config_numerics%getFloat('rtol_crystallitestress',defaultVal=1.0e-6_pReal) num%aTol_crystalliteStress = config_numerics%getFloat('atol_crystallitestress',defaultVal=1.0e-8_pReal) - + num%iJacoLpresiduum = config_numerics%getInt ('ijacolpresiduum', defaultVal=1) - + num%nState = config_numerics%getInt ('nstate', defaultVal=20) num%nStress = config_numerics%getInt ('nstress', defaultVal=40) - + if(num%subStepMinCryst <= 0.0_pReal) call IO_error(301,ext_msg='subStepMinCryst') if(num%subStepSizeCryst <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeCryst') if(num%stepIncreaseCryst <= 0.0_pReal) call IO_error(301,ext_msg='stepIncreaseCryst') @@ -199,12 +199,12 @@ subroutine crystallite_init if(num%rTol_crystalliteState <= 0.0_pReal) call IO_error(301,ext_msg='rTol_crystalliteState') if(num%rTol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='rTol_crystalliteStress') if(num%aTol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='aTol_crystalliteStress') - + if(num%iJacoLpresiduum < 1) call IO_error(301,ext_msg='iJacoLpresiduum') - + if(num%nState < 1) call IO_error(301,ext_msg='nState') if(num%nStress< 1) call IO_error(301,ext_msg='nStress') - + select case(numerics_integrator) case(1) integrateState => integrateStateFPI @@ -217,11 +217,11 @@ subroutine crystallite_init case(5) integrateState => integrateStateRKCK45 end select - + allocate(output_constituent(size(config_phase))) do c = 1, size(config_phase) #if defined(__GFORTRAN__) - allocate(output_constituent(c)%label(1)) + allocate(output_constituent(c)%label(1)) output_constituent(c)%label(1)= 'GfortranBug86277' output_constituent(c)%label = config_phase(c)%getStrings('(output)',defaultVal=output_constituent(c)%label ) if (output_constituent(c)%label (1) == 'GfortranBug86277') output_constituent(c)%label = [character(len=pStringLen)::] @@ -250,16 +250,16 @@ subroutine crystallite_init enddo; enddo enddo !$OMP END PARALLEL DO - + if(any(.not. crystallite_localPlasticity) .and. .not. usePingPong) call IO_error(601) ! exit if nonlocal but no ping-pong ToDo: Why not check earlier? or in nonlocal? - + crystallite_partionedFp0 = crystallite_Fp0 crystallite_partionedFi0 = crystallite_Fi0 crystallite_partionedF0 = crystallite_F0 crystallite_partionedF = crystallite_F0 - + call crystallite_orientations() - + !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) @@ -271,7 +271,7 @@ subroutine crystallite_init enddo enddo !$OMP END PARALLEL DO - + devNull = crystallite_stress() call crystallite_stressTangent @@ -283,7 +283,7 @@ subroutine crystallite_init write(6,'(a42,1x,i10)') ' # of nonlocal constituents: ',count(.not. crystallite_localPlasticity) flush(6) endif - + call debug_info call debug_reset #endif @@ -295,7 +295,7 @@ end subroutine crystallite_init !> @brief calculate stress (P) !-------------------------------------------------------------------------------------------------- function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) - + logical, dimension(discretization_nIP,discretization_nElem) :: crystallite_stress real(pReal), intent(in), optional :: & dummyArgumentToPreventInternalCompilerErrorWithGCC @@ -308,7 +308,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) e, & !< counter in element loop startIP, endIP, & s - + #ifdef DEBUG if (iand(debug_level(debug_crystallite),debug_levelSelective) /= 0 & .and. FEsolving_execElem(1) <= debug_e & @@ -494,7 +494,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) crystallite_stress = .false. elementLooping5: do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) - crystallite_stress(i,e) = all(crystallite_converged(:,i,e)) + crystallite_stress(i,e) = all(crystallite_converged(:,i,e)) enddo enddo elementLooping5 @@ -675,12 +675,12 @@ end subroutine crystallite_stressTangent !> @brief calculates orientations !-------------------------------------------------------------------------------------------------- subroutine crystallite_orientations - + integer & c, & !< counter in integration point component loop i, & !< counter in integration point loop e !< counter in element loop - + !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) @@ -688,7 +688,7 @@ subroutine crystallite_orientations call crystallite_orientation(c,i,e)%fromMatrix(transpose(math_rotationalPart33(crystallite_Fe(1:3,1:3,c,i,e)))) enddo; enddo; enddo !$OMP END PARALLEL DO - + nonlocalPresent: if (any(plasticState%nonLocal)) then !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) @@ -706,7 +706,7 @@ end subroutine crystallite_orientations !> @brief Map 2nd order tensor to reference config !-------------------------------------------------------------------------------------------------- function crystallite_push33ToRef(ipc,ip,el, tensor33) - + real(pReal), dimension(3,3) :: crystallite_push33ToRef real(pReal), dimension(3,3), intent(in) :: tensor33 real(pReal), dimension(3,3) :: T @@ -714,7 +714,7 @@ function crystallite_push33ToRef(ipc,ip,el, tensor33) el, & ip, & ipc - + T = matmul(material_orientation0(ipc,ip,el)%asMatrix(), & ! ToDo: initial orientation correct? transpose(math_inv33(crystallite_subF(1:3,1:3,ipc,ip,el)))) crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T)) @@ -731,11 +731,11 @@ subroutine crystallite_results real(pReal), allocatable, dimension(:,:,:) :: selected_tensors type(rotation), allocatable, dimension(:) :: selected_rotations character(len=256) :: group,lattice_label - + do p=1,size(config_name_phase) group = trim('current/constituent')//'/'//trim(config_name_phase(p))//'/generic' - - call results_closeGroup(results_addGroup(group)) + + call results_closeGroup(results_addGroup(group)) do o = 1, size(output_constituent(p)%label) select case (output_constituent(p)%label(o)) @@ -792,19 +792,19 @@ subroutine crystallite_results end select enddo enddo - + contains !------------------------------------------------------------------------------------------------ !> @brief select tensors for output !------------------------------------------------------------------------------------------------ function select_tensors(dataset,instance) - + integer, intent(in) :: instance real(pReal), dimension(:,:,:,:,:), intent(in) :: dataset real(pReal), allocatable, dimension(:,:,:) :: select_tensors integer :: e,i,c,j - + allocate(select_tensors(3,3,count(material_phaseAt==instance)*homogenization_maxNgrains*discretization_nIP)) j=0 @@ -818,20 +818,20 @@ subroutine crystallite_results enddo enddo enddo - + end function select_tensors - - + + !-------------------------------------------------------------------------------------------------- !> @brief select rotations for output -!-------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- function select_rotations(dataset,instance) - + integer, intent(in) :: instance type(rotation), dimension(:,:,:), intent(in) :: dataset type(rotation), allocatable, dimension(:) :: select_rotations integer :: e,i,c,j - + allocate(select_rotations(count(material_phaseAt==instance)*homogenization_maxNgrains*discretization_nIP)) j=0 @@ -845,7 +845,7 @@ subroutine crystallite_results enddo enddo enddo - + end function select_rotations end subroutine crystallite_results @@ -856,12 +856,12 @@ end subroutine crystallite_results !> intermediate acceleration of the Newton-Raphson correction !-------------------------------------------------------------------------------------------------- logical function integrateStress(ipc,ip,el,timeFraction) - + integer, intent(in):: el, & ! element index ip, & ! integration point index ipc ! grain index real(pReal), optional, intent(in) :: timeFraction ! fraction of timestep - + real(pReal), dimension(3,3):: Fg_new, & ! deformation gradient at end of timestep Fp_new, & ! plastic deformation gradient at end of timestep Fe_new, & ! elastic deformation gradient at end of timestep @@ -916,7 +916,7 @@ logical function integrateStress(ipc,ip,el,timeFraction) jacoCounterLi ! counters to check for Jacobian update external :: & dgesv - + !* be pessimistic integrateStress = .false. #ifdef DEBUG @@ -938,7 +938,7 @@ logical function integrateStress(ipc,ip,el,timeFraction) Lpguess = crystallite_Lp(1:3,1:3,ipc,ip,el) ! take as first guess Liguess = crystallite_Li(1:3,1:3,ipc,ip,el) ! take as first guess Liguess_old = Liguess - + invFp_current = math_inv33(crystallite_subFp0(1:3,1:3,ipc,ip,el)) failedInversionFp: if (all(dEq0(invFp_current))) then #ifdef DEBUG @@ -964,13 +964,13 @@ logical function integrateStress(ipc,ip,el,timeFraction) #endif return endif failedInversionFi - + !* start Li loop with normal step length NiterationStressLi = 0 jacoCounterLi = 0 steplengthLi = 1.0_pReal residuumLi_old = 0.0_pReal - + LiLoop: do NiterationStressLi = NiterationStressLi + 1 LiLoopLimit: if (NiterationStressLi > num%nStress) then @@ -981,18 +981,18 @@ logical function integrateStress(ipc,ip,el,timeFraction) #endif return endif LiLoopLimit - + invFi_new = matmul(invFi_current,math_I3 - dt*Liguess) Fi_new = math_inv33(invFi_new) detInvFi = math_det33(invFi_new) - + !* start Lp loop with normal step length NiterationStressLp = 0 jacoCounterLp = 0 steplengthLp = 1.0_pReal residuumLp_old = 0.0_pReal Lpguess_old = Lpguess - + LpLoop: do NiterationStressLp = NiterationStressLp + 1 LpLoopLimit: if (NiterationStressLp > num%nStress) then @@ -1003,18 +1003,18 @@ logical function integrateStress(ipc,ip,el,timeFraction) #endif return endif LpLoopLimit - + !* calculate (elastic) 2nd Piola--Kirchhoff stress tensor and its tangent from constitutive law - + B = math_I3 - dt*Lpguess Fe = matmul(matmul(A,B), invFi_new) call constitutive_SandItsTangents(S, dS_dFe, dS_dFi, & Fe, Fi_new, ipc, ip, el) ! call constitutive law to calculate 2nd Piola-Kirchhoff stress and its derivative in unloaded configuration - + !* calculate plastic velocity gradient and its tangent from constitutive law call constitutive_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, & S, Fi_new, ipc, ip, el) - + #ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 & .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & @@ -1032,7 +1032,7 @@ logical function integrateStress(ipc,ip,el,timeFraction) aTolLp = max(num%rTol_crystalliteStress * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error num%aTol_crystalliteStress) ! minimum lower cutoff residuumLp = Lpguess - Lp_constitutive - + if (any(IEEE_is_NaN(residuumLp))) then #ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0) & @@ -1160,7 +1160,7 @@ logical function integrateStress(ipc,ip,el,timeFraction) #endif cycle LiLoop endif - + !* calculate Jacobian for correction term if (mod(jacoCounterLi, num%iJacoLpresiduum) == 0) then temp_33 = matmul(matmul(A,B),invFi_current) @@ -1197,11 +1197,11 @@ logical function integrateStress(ipc,ip,el,timeFraction) #endif return endif - + deltaLi = - math_9to33(work) endif jacoCounterLi = jacoCounterLi + 1 - + Liguess = Liguess + steplengthLi * deltaLi #ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 & @@ -1211,7 +1211,7 @@ logical function integrateStress(ipc,ip,el,timeFraction) endif #endif enddo LiLoop - + !* calculate new plastic and elastic deformation gradient invFp_new = matmul(invFp_current,B) invFp_new = invFp_new / math_det33(invFp_new)**(1.0_pReal/3.0_pReal) ! regularize @@ -1302,7 +1302,7 @@ subroutine integrateStateFPI #endif ! store previousDotState and previousDotState2 - + !$OMP PARALLEL DO PRIVATE(p,c) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) @@ -1329,7 +1329,7 @@ subroutine integrateStateFPI call update_dependentState call update_stress(1.0_pReal) call update_dotState(1.0_pReal) - + !$OMP PARALLEL !$OMP DO PRIVATE(sizeDotState,residuum_plastic,residuum_source,zeta,p,c) do e = FEsolving_execElem(1),FEsolving_execElem(2) @@ -1342,7 +1342,7 @@ subroutine integrateStateFPI zeta = damper(plasticState(p)%dotState (:,c), & plasticState(p)%previousDotState (:,c), & plasticState(p)%previousDotState2(:,c)) - + residuum_plastic(1:SizeDotState) = plasticState(p)%state (1:sizeDotState,c) & - plasticState(p)%subState0(1:sizeDotState,c) & - ( plasticState(p)%dotState (:,c) * zeta & @@ -1350,18 +1350,18 @@ subroutine integrateStateFPI ) * crystallite_subdt(g,i,e) plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%state(1:sizeDotState,c) & - - residuum_plastic(1:sizeDotState) + - residuum_plastic(1:sizeDotState) plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * zeta & + plasticState(p)%previousDotState(:,c) * (1.0_pReal - zeta) - + crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState), & plasticState(p)%state(1:sizeDotState,c), & plasticState(p)%aTolState(1:sizeDotState)) - + do s = 1, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - + zeta = damper(sourceState(p)%p(s)%dotState (:,c), & sourceState(p)%p(s)%previousDotState (:,c), & sourceState(p)%p(s)%previousDotState2(:,c)) @@ -1432,12 +1432,12 @@ subroutine integrateStateFPI !> @brief calculate the damping for correction of state and dot state !-------------------------------------------------------------------------------------------------- real(pReal) pure function damper(current,previous,previous2) - + real(pReal), dimension(:), intent(in) ::& current, previous, previous2 - + real(pReal) :: dot_prod12, dot_prod22 - + dot_prod12 = dot_product(current - previous, previous - previous2) dot_prod22 = dot_product(previous - previous2, previous - previous2) if ((dot_product(current,previous) < 0.0_pReal .or. dot_prod12 < 0.0_pReal) .and. dot_prod22 > 0.0_pReal) then @@ -1445,7 +1445,7 @@ subroutine integrateStateFPI else damper = 1.0_pReal endif - + end function damper end subroutine integrateStateFPI @@ -1480,7 +1480,7 @@ subroutine integrateStateAdaptiveEuler c, & s, & sizeDotState - + ! ToDo: MD: once all constitutives use allocate state, attach residuum arrays to the state in case of adaptive Euler real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & @@ -1501,14 +1501,14 @@ subroutine integrateStateAdaptiveEuler if (crystallite_todo(g,i,e)) then p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) sizeDotState = plasticState(p)%sizeDotState - + residuum_plastic(1:sizeDotState,g,i,e) = plasticState(p)%dotstate(1:sizeDotState,c) & * (- 0.5_pReal * crystallite_subdt(g,i,e)) plasticState(p)%state(1:sizeDotState,c) = & plasticState(p)%state(1:sizeDotState,c) + plasticState(p)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) !ToDo: state, partitioned state? do s = 1, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - + residuum_source(1:sizeDotState,s,g,i,e) = sourceState(p)%p(s)%dotstate(1:sizeDotState,c) & * (- 0.5_pReal * crystallite_subdt(g,i,e)) sourceState(p)%p(s)%state(1:sizeDotState,c) = & @@ -1530,17 +1530,17 @@ subroutine integrateStateAdaptiveEuler if (crystallite_todo(g,i,e)) then p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) sizeDotState = plasticState(p)%sizeDotState - + residuum_plastic(1:sizeDotState,g,i,e) = residuum_plastic(1:sizeDotState,g,i,e) & + 0.5_pReal * plasticState(p)%dotState(:,c) * crystallite_subdt(g,i,e) - + crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState,g,i,e), & plasticState(p)%state(1:sizeDotState,c), & plasticState(p)%aTolState(1:sizeDotState)) do s = 1, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - + residuum_source(1:sizeDotState,s,g,i,e) = & residuum_source(1:sizeDotState,s,g,i,e) + 0.5_pReal * sourceState(p)%p(s)%dotState(:,c) * crystallite_subdt(g,i,e) @@ -1549,13 +1549,13 @@ subroutine integrateStateAdaptiveEuler sourceState(p)%p(s)%state(1:sizeDotState,c), & sourceState(p)%p(s)%aTolState(1:sizeDotState)) enddo - + endif enddo; enddo; enddo !$OMP END PARALLEL DO if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck - + end subroutine integrateStateAdaptiveEuler @@ -1720,23 +1720,23 @@ subroutine integrateStateRKCK45 do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) if (crystallite_todo(g,i,e)) then p = material_phaseAt(g,e); cc = material_phaseMemberAt(g,i,e) - + sizeDotState = plasticState(p)%sizeDotState - + plasticState(p)%RKCK45dotState(6,:,cc) = plasticState (p)%dotState(:,cc) - + residuum_plastic(1:sizeDotState,g,i,e) = & matmul(transpose(plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,cc)),DB) & ! why transpose? Better to transpose constant DB * crystallite_subdt(g,i,e) - + plasticState(p)%dotState(:,cc) = & matmul(transpose(plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,cc)), B) ! why transpose? Better to transpose constant B - + do s = 1, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - + sourceState(p)%p(s)%RKCK45dotState(6,:,cc) = sourceState(p)%p(s)%dotState(:,cc) - + residuum_source(1:sizeDotState,s,g,i,e) = & matmul(transpose(sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,cc)),DB) & * crystallite_subdt(g,i,e) @@ -1744,13 +1744,13 @@ subroutine integrateStateRKCK45 sourceState(p)%p(s)%dotState(:,cc) = & matmul(transpose(sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,cc)),B) enddo - + endif enddo; enddo; enddo !$OMP END PARALLEL DO call update_state(1.0_pReal) - + ! --- relative residui and state convergence --- !$OMP PARALLEL DO PRIVATE(sizeDotState,p,cc) @@ -1759,16 +1759,16 @@ subroutine integrateStateRKCK45 do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) if (crystallite_todo(g,i,e)) then p = material_phaseAt(g,e); cc = material_phaseMemberAt(g,i,e) - + sizeDotState = plasticState(p)%sizeDotState - + crystallite_todo(g,i,e) = converged(residuum_plastic(1:sizeDotState,g,i,e), & plasticState(p)%state(1:sizeDotState,cc), & plasticState(p)%aTolState(1:sizeDotState)) do s = 1, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - + crystallite_todo(g,i,e) = & crystallite_todo(g,i,e) .and. converged(residuum_source(1:sizeDotState,s,g,i,e), & sourceState(p)%p(s)%state(1:sizeDotState,cc), & @@ -1783,7 +1783,7 @@ subroutine integrateStateRKCK45 call update_stress(1.0_pReal) call setConvergenceFlag if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck - + end subroutine integrateStateRKCK45 @@ -1792,7 +1792,7 @@ end subroutine integrateStateRKCK45 !> @detail one non-converged nonlocal sets all other nonlocals to non-converged to trigger cut back !-------------------------------------------------------------------------------------------------- subroutine nonlocalConvergenceCheck - + if (any(.not. crystallite_converged .and. .not. crystallite_localPlasticity)) & ! any non-local not yet converged (or broken)... where( .not. crystallite_localPlasticity) crystallite_converged = .false. @@ -1810,7 +1810,7 @@ subroutine setConvergenceFlag e, & !< element index in element loop i, & !< integration point index in ip loop g !< grain index in grain loop - + !OMP DO PARALLEL PRIVATE do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) @@ -1826,7 +1826,7 @@ end subroutine setConvergenceFlag !> @brief determines whether a point is converged !-------------------------------------------------------------------------------------------------- logical pure function converged(residuum,state,aTol) - + real(pReal), intent(in), dimension(:) ::& residuum, state, aTol real(pReal) :: & @@ -1949,11 +1949,11 @@ subroutine update_dotState(timeFraction) g, & !< grain index in grain loop p, & c, & - s + s logical :: & NaN, & nonlocalStop - + nonlocalStop = .false. !$OMP PARALLEL DO PRIVATE (p,c,NaN) @@ -1963,9 +1963,9 @@ subroutine update_dotState(timeFraction) !$OMP FLUSH(nonlocalStop) if ((crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) .and. .not. nonlocalStop) then call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_Fe, & + crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & - crystallite_Fp, & + crystallite_partionedFp0, & crystallite_subdt(g,i,e)*timeFraction, g,i,e) p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c))) @@ -1980,7 +1980,7 @@ subroutine update_dotState(timeFraction) enddo; enddo; enddo !$OMP END PARALLEL DO - if (nonlocalStop) crystallite_todo = crystallite_todo .and. crystallite_localPlasticity + if (nonlocalStop) crystallite_todo = crystallite_todo .and. crystallite_localPlasticity end subroutine update_DotState @@ -1995,11 +1995,11 @@ subroutine update_deltaState mySize, & myOffset, & c, & - s + s logical :: & NaN, & nonlocalStop - + nonlocalStop = .false. !$OMP PARALLEL DO PRIVATE(p,c,myOffset,mySize,NaN) @@ -2016,23 +2016,23 @@ subroutine update_deltaState myOffset = plasticState(p)%offsetDeltaState mySize = plasticState(p)%sizeDeltaState NaN = any(IEEE_is_NaN(plasticState(p)%deltaState(1:mySize,c))) - + if (.not. NaN) then - + plasticState(p)%state(myOffset + 1: myOffset + mySize,c) = & plasticState(p)%state(myOffset + 1: myOffset + mySize,c) + plasticState(p)%deltaState(1:mySize,c) do s = 1, phase_Nsources(p) myOffset = sourceState(p)%p(s)%offsetDeltaState mySize = sourceState(p)%p(s)%sizeDeltaState NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(s)%deltaState(1:mySize,c))) - + if (.not. NaN) then sourceState(p)%p(s)%state(myOffset + 1:myOffset + mySize,c) = & sourceState(p)%p(s)%state(myOffset + 1:myOffset + mySize,c) + sourceState(p)%p(s)%deltaState(1:mySize,c) endif enddo endif - + crystallite_todo(g,i,e) = .not. NaN if (.not. crystallite_todo(g,i,e)) then ! if state jump fails, then convergence is broken crystallite_converged(g,i,e) = .false. @@ -2042,7 +2042,7 @@ subroutine update_deltaState enddo; enddo; enddo !$OMP END PARALLEL DO if (nonlocalStop) crystallite_todo = crystallite_todo .and. crystallite_localPlasticity - + end subroutine update_deltaState From 5d4d1dcf9a175500e5ba2207ef20f6e1e39d5044 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 7 Feb 2020 12:41:01 +0100 Subject: [PATCH 4/7] all nonlocal parts are fully explicit i.e. they are based on converged (partioned0) states --- src/constitutive.f90 | 10 +++++----- src/constitutive_plastic_nonlocal.f90 | 27 ++++++++++++++------------- src/crystallite.f90 | 4 ++-- 3 files changed, 21 insertions(+), 20 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index f89c08be8..6f3c91fff 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -232,12 +232,12 @@ module constitutive of end subroutine plastic_disloUCLA_dependentState - module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) + module subroutine plastic_nonlocal_dependentState(F, Fp, ip, el) integer, intent(in) :: & ip, & el real(pReal), dimension(3,3), intent(in) :: & - Fe, & + F, & Fp end subroutine plastic_nonlocal_dependentState @@ -412,14 +412,14 @@ end function constitutive_homogenizedC !-------------------------------------------------------------------------------------------------- !> @brief calls microstructure function of the different constitutive models !-------------------------------------------------------------------------------------------------- -subroutine constitutive_microstructure(Fe, Fp, ipc, ip, el) +subroutine constitutive_microstructure(F, Fp, ipc, ip, el) integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element real(pReal), intent(in), dimension(3,3) :: & - Fe, & !< elastic deformation gradient + F, & !< elastic deformation gradient Fp !< plastic deformation gradient integer :: & ho, & !< homogenization @@ -439,7 +439,7 @@ subroutine constitutive_microstructure(Fe, Fp, ipc, ip, el) instance = phase_plasticityInstance(material_phaseAt(ipc,el)) call plastic_disloUCLA_dependentState(instance,of) case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_dependentState (Fe,Fp,ip,el) + call plastic_nonlocal_dependentState (F,Fp,ip,el) end select plasticityType end subroutine constitutive_microstructure diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index 0a51e78d9..8bf6ac6fc 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -716,13 +716,13 @@ end subroutine plastic_nonlocal_init !-------------------------------------------------------------------------------------------------- !> @brief calculates quantities characterizing the microstructure !-------------------------------------------------------------------------------------------------- -module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) +module subroutine plastic_nonlocal_dependentState(F, Fp, ip, el) integer, intent(in) :: & ip, & el real(pReal), dimension(3,3), intent(in) :: & - Fe, & + F, & Fp integer :: & @@ -765,7 +765,8 @@ module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) rho_scr_delta real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),10) :: & rho, & - rho_neighbor + rho0, & + rho_neighbor0 real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))), & totalNslip(phase_plasticityInstance(material_phaseAt(1,el)))) :: & myInteractionMatrix ! corrected slip interaction matrix @@ -795,7 +796,7 @@ module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) ! coefficients are corrected for the line tension effect ! (see Kubin,Devincre,Hoc; 2008; Modeling dislocation storage rates and mean free paths in face-centered cubic crystals) - if (lattice_structure(ph) == LATTICE_bcc_ID .or. lattice_structure(ph) == LATTICE_fcc_ID) then ! only fcc and bcc + if (lattice_structure(ph) == LATTICE_bcc_ID .or. lattice_structure(ph) == LATTICE_fcc_ID) then do s = 1,ns correction = ( 1.0_pReal - prm%linetensionEffect & + prm%linetensionEffect & @@ -819,13 +820,13 @@ module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) ! ToDo: MD: this is most likely only correct for F_i = I !################################################################################################# - + rho0 = getRho0(instance,of,ip,el) if (.not. phase_localPlasticity(ph) .and. prm%shortRangeStressCorrection) then - invFe = math_inv33(Fe) invFp = math_inv33(Fp) + invFe = matmul(Fp,math_inv33(F)) - rho_edg_delta = rho(:,mob_edg_pos) - rho(:,mob_edg_neg) - rho_scr_delta = rho(:,mob_scr_pos) - rho(:,mob_scr_neg) + rho_edg_delta = rho0(:,mob_edg_pos) - rho0(:,mob_edg_neg) + rho_scr_delta = rho0(:,mob_scr_pos) - rho0(:,mob_scr_neg) rhoExcess(1,1:ns) = rho_edg_delta rhoExcess(2,1:ns) = rho_scr_delta @@ -845,13 +846,13 @@ module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) if (neighbor_instance == instance) then nRealNeighbors = nRealNeighbors + 1.0_pReal - rho_neighbor = getRho(instance,no,neighbor_ip,neighbor_el) + rho_neighbor0 = getRho0(instance,no,neighbor_ip,neighbor_el) - rho_edg_delta_neighbor(:,n) = rho_neighbor(:,mob_edg_pos) - rho_neighbor(:,mob_edg_neg) - rho_scr_delta_neighbor(:,n) = rho_neighbor(:,mob_scr_pos) - rho_neighbor(:,mob_scr_neg) + rho_edg_delta_neighbor(:,n) = rho_neighbor0(:,mob_edg_pos) - rho_neighbor0(:,mob_edg_neg) + rho_scr_delta_neighbor(:,n) = rho_neighbor0(:,mob_scr_pos) - rho_neighbor0(:,mob_scr_neg) - neighbor_rhoTotal(1,:,n) = sum(abs(rho_neighbor(:,edg)),2) - neighbor_rhoTotal(2,:,n) = sum(abs(rho_neighbor(:,scr)),2) + neighbor_rhoTotal(1,:,n) = sum(abs(rho_neighbor0(:,edg)),2) + neighbor_rhoTotal(2,:,n) = sum(abs(rho_neighbor0(:,scr)),2) connection_latticeConf(1:3,n) = matmul(invFe, discretization_IPcoords(1:3,neighbor_el+neighbor_ip-1) & - discretization_IPcoords(1:3,el+neighbor_ip-1)) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 53164c8ee..83ef3fefc 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -264,8 +264,8 @@ subroutine crystallite_init do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do c = 1,homogenization_Ngrains(material_homogenizationAt(e)) - call constitutive_microstructure(crystallite_Fe(1:3,1:3,c,i,e), & - crystallite_Fp(1:3,1:3,c,i,e), & + call constitutive_microstructure(crystallite_partionedF0(1:3,1:3,c,i,e), & + crystallite_partionedFp0(1:3,1:3,c,i,e), & c,i,e) ! update dependent state variables to be consistent with basic states enddo enddo From 9c2f9bfceecb7072db60e93f5633c3b4d78a6e57 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 7 Feb 2020 12:44:35 +0100 Subject: [PATCH 5/7] use relaxed test --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index 64432754c..9afd3d55e 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 64432754ce3c590c882cf4987695539cee524ee8 +Subproject commit 9afd3d55e0ed90fd1a001180d15c4999eb8b4d4a From a306e473efb358bb169c8ff53c243efc6f1c2465 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 11 Feb 2020 05:41:10 +0100 Subject: [PATCH 6/7] use rhoSgl0 (converged situation) --- src/constitutive_plastic_nonlocal.f90 | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index 8bf6ac6fc..13d629e2b 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -1363,6 +1363,7 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature, & s !< index of my current slip system real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),10) :: & rho, & + rho0, & !< dislocation density at beginning of time step rhoDot, & !< density evolution rhoDotMultiplication, & !< density evolution by multiplication rhoDotFlux, & !< density evolution by flux @@ -1371,8 +1372,8 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature, & rhoDotThermalAnnihilation !< density evolution by thermal annihilation real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),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_rhoSgl !< single dislocation densities of central ip (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(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),4) :: & v, & !< current dislocation glide velocity v0, & @@ -1428,9 +1429,11 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature, & tau = 0.0_pReal gdot = 0.0_pReal - rho = getRho(instance,o,ip,el) + rho = getRho(instance,o,ip,el) rhoSgl = rho(:,sgl) rhoDip = rho(:,dip) + rho0 = getRho0(instance,o,ip,el) + my_rhoSgl0 = rho0(:,sgl) forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(p)%state(iV (s,t,instance),o) @@ -1626,8 +1629,6 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature, & endif leavingFlux: if (considerLeavingFlux) then - my_rhoSgl = rhoSgl - normal_me2neighbor_defConf = math_det33(Favg) & * matmul(math_inv33(transpose(Favg)), & IPareaNormal(1:3,n,ip,el)) ! calculate the normal of the interface in (average) deformed configuration (pointing from me to my neighbor!!!) @@ -1644,7 +1645,7 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature, & 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_rhoSgl(s,t) * v0(s,t) & + 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) & From 61810e6a08cc5c1dd362931eb161bcb60749c820 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 11 Feb 2020 05:54:19 +0100 Subject: [PATCH 7/7] use tighter tolerance (works again) --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index 9afd3d55e..9dc7065be 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 9afd3d55e0ed90fd1a001180d15c4999eb8b4d4a +Subproject commit 9dc7065beedf2a097aa60a656bbfa52e55d7147c