!-------------------------------------------------------------------------------------------------- !> @author Christoph Kords, Max-Planck-Institut für Eisenforschung GmbH !> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @brief material subroutine for plasticity including dislocation flux !-------------------------------------------------------------------------------------------------- module plastic_nonlocal use prec, only: & pReal, & pInt implicit none private real(pReal), parameter, private :: & KB = 1.38e-23_pReal !< Physical parameter, Boltzmann constant in J/Kelvin integer(pInt), dimension(:,:), allocatable, target, public :: & plastic_nonlocal_sizePostResult !< size of each post result output character(len=64), dimension(:,:), allocatable, target, public :: & plastic_nonlocal_output !< name of each post result output integer, dimension(8), parameter :: & sgl = [1,2,3,4,5,6,7,8] integer, dimension(2), parameter :: & dip = [9,10] integer, dimension(5), parameter :: & edg = [1,2,5,6,9], & scr = [3,4,7,8,10] integer, dimension(4), parameter :: & mob = [1,2,3,4], & imm = [5,6,7,8], & pos = sgl(1:7:2), & neg = sgl(2:8:2), & sgl_edg = edg(1:4), & sgl_scr = scr(1:4) integer, parameter :: & mob_edg_pos = 1, & mob_edg_neg = 2, & mob_scr_pos = 3, & mob_scr_neg = 4 integer(pInt), dimension(:,:), allocatable, private :: & iRhoF !< state indices for forest density integer(pInt), dimension(:,:,:), allocatable, private :: & iRhoU, & !< state indices for unblocked density iRhoB, & !< state indices for blocked density iRhoD, & !< state indices for dipole density iV, & !< state indices for dislcation velocities iD !< state indices for stable dipole height integer(pInt), dimension(:), allocatable, public, protected :: & totalNslip !< total number of active slip systems for each instance real(pReal), dimension(:,:,:,:,:,:), allocatable, private :: & compatibility !< slip system compatibility between me and my neighbors enum, bind(c) enumerator :: undefined_ID, & rho_sgl_mob_edg_pos_ID, & rho_sgl_mob_edg_neg_ID, & rho_sgl_mob_scr_pos_ID, & rho_sgl_mob_scr_neg_ID, & rho_sgl_imm_edg_pos_ID, & rho_sgl_imm_edg_neg_ID, & rho_sgl_imm_scr_pos_ID, & rho_sgl_imm_scr_neg_ID, & rho_dip_edg_ID, & rho_dip_scr_ID, & rho_forest_ID, & shearrate_ID, & resolvedstress_ID, & resolvedstress_external_ID, & resolvedstress_back_ID, & resistance_ID, & rho_dot_sgl_ID, & rho_dot_sgl_mobile_ID, & rho_dot_dip_ID, & rho_dot_gen_ID, & rho_dot_gen_edge_ID, & rho_dot_gen_screw_ID, & rho_dot_sgl2dip_edge_ID, & rho_dot_sgl2dip_screw_ID, & rho_dot_ann_ath_ID, & rho_dot_ann_the_edge_ID, & rho_dot_ann_the_screw_ID, & rho_dot_edgejogs_ID, & rho_dot_flux_mobile_ID, & rho_dot_flux_edge_ID, & rho_dot_flux_screw_ID, & velocity_edge_pos_ID, & velocity_edge_neg_ID, & velocity_screw_pos_ID, & velocity_screw_neg_ID, & maximumdipoleheight_edge_ID, & maximumdipoleheight_screw_ID, & accumulatedshear_ID end enum type, private :: tParameters !< container type for internal constitutive parameters real(pReal) :: & atomicVolume, & !< atomic volume Dsd0, & !< prefactor for self-diffusion coefficient selfDiffusionEnergy, & !< activation enthalpy for diffusion aTolRho, & !< absolute tolerance for dislocation density in state integration aTolShear, & !< absolute tolerance for accumulated shear in state integration significantRho, & !< density considered significant significantN, & !< number of dislocations considered significant doublekinkwidth, & !< width of a doubkle kink in multiples of the burgers vector length b solidSolutionEnergy, & !< activation energy for solid solution in J solidSolutionSize, & !< solid solution obstacle size in multiples of the burgers vector length solidSolutionConcentration, & !< concentration of solid solution in atomic parts p, & !< parameter for kinetic law (Kocks,Argon,Ashby) q, & !< parameter for kinetic law (Kocks,Argon,Ashby) viscosity, & !< viscosity for dislocation glide in Pa s fattack, & !< attack frequency in Hz rhoSglScatter, & !< standard deviation of scatter in initial dislocation density surfaceTransmissivity, & !< transmissivity at free surface grainboundaryTransmissivity, & !< transmissivity at grain boundary (identified by different texture) CFLfactor, & !< safety factor for CFL flux condition fEdgeMultiplication, & !< factor that determines how much edge dislocations contribute to multiplication (0...1) rhoSglRandom, & rhoSglRandomBinning, & linetensionEffect, & edgeJogFactor, & mu, & nu real(pReal), dimension(:), allocatable :: & minDipoleHeight_edge, & !< minimum stable edge dipole height minDipoleHeight_screw, & !< minimum stable screw dipole height peierlsstress_edge, & peierlsstress_screw, & rhoSglEdgePos0, & !< initial edge_pos dislocation density rhoSglEdgeNeg0, & !< initial edge_neg dislocation density rhoSglScrewPos0, & !< initial screw_pos dislocation density rhoSglScrewNeg0, & !< initial screw_neg dislocation density rhoDipEdge0, & !< initial edge dipole dislocation density rhoDipScrew0,& !< initial screw dipole dislocation density lambda0, & !< mean free path prefactor for each burgers !< absolute length of burgers vector [m] real(pReal), dimension(:,:), allocatable :: & slip_normal, & slip_direction, & slip_transverse, & minDipoleHeight, & ! edge and screw peierlsstress, & ! edge and screw interactionSlipSlip ,& !< coefficients for slip-slip interaction forestProjection_Edge, & !< matrix of forest projections of edge dislocations forestProjection_Screw !< matrix of forest projections of screw dislocations real(pReal), dimension(:), allocatable, private :: & nonSchmidCoeff integer(pInt) :: totalNslip real(pReal), dimension(:,:,:), allocatable, private :: & Schmid, & !< Schmid contribution nonSchmid_pos, & nonSchmid_neg !< combined projection of Schmid and non-Schmid contributions to the resolved shear stress (only for screws) integer(pInt) , dimension(:) ,allocatable , public:: & Nslip,& colinearSystem !< colinear system to the active slip system (only valid for fcc!) logical, private :: & 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, private :: tNonlocalMicrostructure real(pReal), allocatable, dimension(:,:) :: & tau_Threshold, & tau_Back end type tNonlocalMicrostructure type, private :: tOutput !< container type for storage of output results real(pReal), dimension(:,:), allocatable, private :: & rhoDotEdgeJogs real(pReal), dimension(:,:,:), allocatable, private :: & rhoDotFlux, & rhoDotMultiplication, & rhoDotSingle2DipoleGlide, & rhoDotAthermalAnnihilation, & rhoDotThermalAnnihilation end type type, private :: tNonlocalState real(pReal), pointer, dimension(:,:) :: & rho, & ! < all dislocations rhoSgl, & rhoSglMobile, & ! iRhoU rhoSglEdgeMobile, & rho_sgl_mob_edg_pos, & rho_sgl_mob_edg_neg, & rhoSglScrewMobile, & rho_sgl_mob_scr_pos, & rho_sgl_mob_scr_neg, & rhoSglImmobile, & ! iRhoB rhoSglEdgeImmobile, & rho_sgl_imm_edg_pos, & rho_sgl_imm_edg_neg, & rhoSglScrewImmobile, & rho_sgl_imm_scr_pos, & rho_sgl_imm_scr_neg, & rhoSglPos, & rhoSglMobilePos, & rhoSglImmobilePos, & rhoSglNeg, & rhoSglMobileNeg, & rhoSglImmobileNeg, & rhoDip, & ! iRhoD rho_dip_edg, & rho_dip_scr, & rhoSglScrew, & rhoSglEdge, & rho_forest, & accumulatedshear end type tNonlocalState type(tNonlocalState), allocatable, dimension(:), private :: & deltaState, & dotState, & state type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) type(tOutput), dimension(:), allocatable, private :: results type(tNonlocalMicrostructure), dimension(:), allocatable, private :: microstructure integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: & plastic_nonlocal_outputID !< ID of each post result output public :: & plastic_nonlocal_init, & plastic_nonlocal_dependentState, & plastic_nonlocal_LpAndItsTangent, & plastic_nonlocal_dotState, & plastic_nonlocal_deltaState, & plastic_nonlocal_updateCompatibility, & plastic_nonlocal_postResults, & plastic_nonlocal_results private :: & plastic_nonlocal_kinetics contains !-------------------------------------------------------------------------------------------------- !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine plastic_nonlocal_init use prec, only: & dEq0, dNeq0, dEq use math, only: & math_expand, math_cross use IO, only: & IO_error use debug, only: & debug_level, & debug_constitutive, & debug_levelBasic use mesh, only: & theMesh use material, only: & phase_plasticity, & phase_plasticityInstance, & phase_Noutput, & PLASTICITY_NONLOCAL_label, & PLASTICITY_NONLOCAL_ID, & plasticState, & material_phase, & material_allocatePlasticState use config use lattice implicit none character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] integer(pInt), dimension(0), parameter :: emptyIntArray = [integer(pInt)::] real(pReal), dimension(0), parameter :: emptyRealArray = [real(pReal)::] integer(pInt) :: & maxNinstances, & p, i, & l, & s1, s2, & s, & ! index of my slip system t, & ! index of dislocation type c ! index of dislocation character integer(pInt) :: sizeState, sizeDotState,sizeDependentState, sizeDeltaState integer(kind(undefined_ID)) :: & outputID character(len=512) :: & extmsg = '', & structure character(len=65536), dimension(:), allocatable :: outputs integer(pInt) :: NofMyPhase write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_NONLOCAL_label//' init -+>>>' write(6,'(/,a)') ' Reuber et al., Acta Materialia 71:333–348, 2014' write(6,'(a)') ' https://doi.org/10.1016/j.actamat.2014.03.012' write(6,'(/,a)') ' Kords, Dissertation RWTH Aachen, 2014' write(6,'(a)') ' http://publications.rwth-aachen.de/record/229993' maxNinstances = count(phase_plasticity == PLASTICITY_NONLOCAL_ID) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & write(6,'(a16,1x,i5,/)') '# instances:',maxNinstances allocate(param(maxNinstances)) allocate(state(maxNinstances)) allocate(dotState(maxNinstances)) allocate(deltaState(maxNinstances)) allocate(microstructure(maxNinstances)) allocate(results(maxNinstances)) allocate(plastic_nonlocal_sizePostResult(maxval(phase_Noutput), maxNinstances), source=0_pInt) allocate(plastic_nonlocal_output(maxval(phase_Noutput), maxNinstances)) plastic_nonlocal_output = '' allocate(plastic_nonlocal_outputID(maxval(phase_Noutput), maxNinstances), source=undefined_ID) allocate(totalNslip(maxNinstances), source=0_pInt) do p=1_pInt, size(config_phase) if (phase_plasticity(p) /= PLASTICITY_NONLOCAL_ID) cycle associate(prm => param(phase_plasticityInstance(p)), & dot => dotState(phase_plasticityInstance(p)), & stt => state(phase_plasticityInstance(p)), & del => deltaState(phase_plasticityInstance(p)), & res => results(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 prm%mu = lattice_mu(p) prm%nu = lattice_nu(p) prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) prm%totalNslip = sum(prm%Nslip) slipActive: if (prm%totalNslip > 0_pInt) then prm%Schmid = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),& config%getFloat('c/a',defaultVal=0.0_pReal)) if(trim(config%getString('lattice_structure')) == 'bcc') then prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',& defaultVal = emptyRealArray) prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt) prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1_pInt) else prm%nonSchmid_pos = prm%Schmid prm%nonSchmid_neg = prm%Schmid endif prm%interactionSlipSlip = lattice_interaction_SlipBySlip(prm%Nslip, & config%getFloats('interaction_slipslip'), & config%getString('lattice_structure')) prm%forestProjection_edge = lattice_forestProjection_edge (prm%Nslip,config%getString('lattice_structure'),& config%getFloat('c/a',defaultVal=0.0_pReal)) prm%forestProjection_screw = lattice_forestProjection_screw(prm%Nslip,config%getString('lattice_structure'),& config%getFloat('c/a',defaultVal=0.0_pReal)) prm%slip_direction = lattice_slip_direction (prm%Nslip,config%getString('lattice_structure'),& config%getFloat('c/a',defaultVal=0.0_pReal)) prm%slip_transverse = lattice_slip_transverse(prm%Nslip,config%getString('lattice_structure'),& config%getFloat('c/a',defaultVal=0.0_pReal)) prm%slip_normal = lattice_slip_normal (prm%Nslip,config%getString('lattice_structure'),& config%getFloat('c/a',defaultVal=0.0_pReal)) ! collinear systems (only for octahedral slip systems in fcc) allocate(prm%colinearSystem(prm%totalNslip), source = -1_pInt) do s1 = 1_pInt, prm%totalNslip do s2 = 1_pInt, prm%totalNslip if (all(dEq0 (math_cross(prm%slip_direction(1:3,s1),prm%slip_direction(1:3,s2)))) .and. & any(dNeq0(math_cross(prm%slip_normal (1:3,s1),prm%slip_normal (1:3,s2))))) & 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) prm%minDipoleHeight_screw = math_expand(prm%minDipoleHeight_screw,prm%Nslip) 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) prm%peierlsstress_screw = math_expand(prm%peierlsstress_screw,prm%Nslip) allocate(prm%peierlsstress(prm%totalNslip,2)) prm%peierlsstress(:,1) = prm%peierlsstress_edge prm%peierlsstress(:,2) = prm%peierlsstress_screw 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') prm%linetensionEffect = config%getFloat('linetension') prm%edgeJogFactor = config%getFloat('edgejog')!,'edgejogs' prm%doublekinkwidth = config%getFloat('doublekinkwidth') 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') prm%fattack = config%getFloat('attackfrequency') ! ToDo: discuss logic prm%rhoSglScatter = config%getFloat('rhosglscatter') prm%rhoSglRandom = config%getFloat('rhosglrandom',0.0_pReal) 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%getInt('shortrangestresscorrection',defaultVal=0_pInt ) > 0_pInt ! ToDo: use /flag/ type key !-------------------------------------------------------------------------------------------------- ! 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%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%surfaceTransmissivity < 0.0_pReal .or. prm%surfaceTransmissivity > 1.0_pReal) & extmsg = trim(extmsg)//' surfaceTransmissivity' if (prm%fEdgeMultiplication < 0.0_pReal .or. prm%fEdgeMultiplication > 1.0_pReal) & extmsg = trim(extmsg)//' fEdgeMultiplication' endif slipActive !-------------------------------------------------------------------------------------------------- ! output pararameters outputs = config%getStrings('(output)',defaultVal=emptyStringArray) allocate(prm%outputID(0)) do i=1_pInt, size(outputs) outputID = undefined_ID select case(trim(outputs(i))) case ('rho_sgl_edge_pos_mobile') outputID = merge(rho_sgl_mob_edg_pos_ID,undefined_ID,prm%totalNslip>0_pInt) case ('rho_sgl_edge_neg_mobile') outputID = merge(rho_sgl_mob_edg_neg_ID,undefined_ID,prm%totalNslip>0_pInt) case ('rho_sgl_screw_pos_mobile') outputID = merge(rho_sgl_mob_scr_pos_ID,undefined_ID,prm%totalNslip>0_pInt) case ('rho_sgl_screw_neg_mobile') outputID = merge(rho_sgl_mob_scr_neg_ID,undefined_ID,prm%totalNslip>0_pInt) case ('rho_sgl_edge_pos_immobile') outputID = merge(rho_sgl_imm_edg_pos_ID,undefined_ID,prm%totalNslip>0_pInt) case ('rho_sgl_edge_neg_immobile') outputID = merge(rho_sgl_imm_edg_neg_ID,undefined_ID,prm%totalNslip>0_pInt) case ('rho_sgl_screw_pos_immobile') outputID = merge(rho_sgl_imm_scr_pos_ID,undefined_ID,prm%totalNslip>0_pInt) case ('rho_sgl_screw_neg_immobile') outputID = merge(rho_sgl_imm_scr_neg_ID,undefined_ID,prm%totalNslip>0_pInt) case ('rho_dip_edge') outputID = merge(rho_dip_edg_ID,undefined_ID,prm%totalNslip>0_pInt) case ('rho_dip_screw') outputID = merge(rho_dip_scr_ID,undefined_ID,prm%totalNslip>0_pInt) case ('rho_forest') outputID = merge(rho_forest_ID,undefined_ID,prm%totalNslip>0_pInt) case ('shearrate') outputID = merge(shearrate_ID,undefined_ID,prm%totalNslip>0_pInt) case ('resolvedstress') outputID = merge(resolvedstress_ID,undefined_ID,prm%totalNslip>0_pInt) case ('resolvedstress_external') outputID = merge(resolvedstress_external_ID,undefined_ID,prm%totalNslip>0_pInt) case ('resolvedstress_back') outputID = merge(resolvedstress_back_ID,undefined_ID,prm%totalNslip>0_pInt) case ('resistance') outputID = merge(resistance_ID,undefined_ID,prm%totalNslip>0_pInt) case ('rho_dot_sgl') outputID = merge(rho_dot_sgl_ID,undefined_ID,prm%totalNslip>0_pInt) case ('rho_dot_sgl_mobile') outputID = merge(rho_dot_sgl_mobile_ID,undefined_ID,prm%totalNslip>0_pInt) case ('rho_dot_dip') outputID = merge(rho_dot_dip_ID,undefined_ID,prm%totalNslip>0_pInt) case ('rho_dot_gen') outputID = merge(rho_dot_gen_ID,undefined_ID,prm%totalNslip>0_pInt) case ('rho_dot_gen_edge') outputID = merge(rho_dot_gen_edge_ID,undefined_ID,prm%totalNslip>0_pInt) case ('rho_dot_gen_screw') outputID = merge(rho_dot_gen_screw_ID,undefined_ID,prm%totalNslip>0_pInt) case ('rho_dot_sgl2dip_edge') outputID = merge(rho_dot_sgl2dip_edge_ID,undefined_ID,prm%totalNslip>0_pInt) case ('rho_dot_sgl2dip_screw') outputID = merge(rho_dot_sgl2dip_screw_ID,undefined_ID,prm%totalNslip>0_pInt) case ('rho_dot_ann_ath') outputID = merge(rho_dot_ann_ath_ID,undefined_ID,prm%totalNslip>0_pInt) case ('rho_dot_ann_the_edge') outputID = merge(rho_dot_ann_the_edge_ID,undefined_ID,prm%totalNslip>0_pInt) case ('rho_dot_ann_the_screw') outputID = merge(rho_dot_ann_the_screw_ID,undefined_ID,prm%totalNslip>0_pInt) case ('rho_dot_edgejogs') outputID = merge(rho_dot_edgejogs_ID,undefined_ID,prm%totalNslip>0_pInt) case ('rho_dot_flux_mobile') outputID = merge(rho_dot_flux_mobile_ID,undefined_ID,prm%totalNslip>0_pInt) case ('rho_dot_flux_edge') outputID = merge(rho_dot_flux_edge_ID,undefined_ID,prm%totalNslip>0_pInt) case ('rho_dot_flux_screw') outputID = merge(rho_dot_flux_screw_ID,undefined_ID,prm%totalNslip>0_pInt) case ('velocity_edge_pos') outputID = merge(velocity_edge_pos_ID,undefined_ID,prm%totalNslip>0_pInt) case ('velocity_edge_neg') outputID = merge(velocity_edge_neg_ID,undefined_ID,prm%totalNslip>0_pInt) case ('velocity_screw_pos') outputID = merge(velocity_screw_pos_ID,undefined_ID,prm%totalNslip>0_pInt) case ('velocity_screw_neg') outputID = merge(velocity_screw_neg_ID,undefined_ID,prm%totalNslip>0_pInt) case ('maximumdipoleheight_edge') outputID = merge(maximumdipoleheight_edge_ID,undefined_ID,prm%totalNslip>0_pInt) case ('maximumdipoleheight_screw') outputID = merge(maximumdipoleheight_screw_ID,undefined_ID,prm%totalNslip>0_pInt) case ('accumulatedshear','accumulated_shear') outputID = merge(accumulatedshear_ID,undefined_ID,prm%totalNslip>0_pInt) end select if (outputID /= undefined_ID) then plastic_nonlocal_output(i,phase_plasticityInstance(p)) = outputs(i) plastic_nonlocal_sizePostResult(i,phase_plasticityInstance(p)) = prm%totalNslip prm%outputID = [prm%outputID , outputID] endif enddo !-------------------------------------------------------------------------------------------------- ! allocate state arrays NofMyPhase=count(material_phase==p) sizeDotState = int(size([ 'rhoSglEdgePosMobile ','rhoSglEdgeNegMobile ', & 'rhoSglScrewPosMobile ','rhoSglScrewNegMobile ', & 'rhoSglEdgePosImmobile ','rhoSglEdgeNegImmobile ', & 'rhoSglScrewPosImmobile','rhoSglScrewNegImmobile', & 'rhoDipEdge ','rhoDipScrew ', & 'accumulatedshear ' ]),pInt) * prm%totalNslip !< "basic" microstructural state variables that are independent from other state variables sizeDependentState = int(size([ 'rhoForest ']),pInt) * prm%totalNslip !< microstructural state variables that depend on other state variables sizeState = sizeDotState + sizeDependentState & + int(size([ 'velocityEdgePos ','velocityEdgeNeg ', & 'velocityScrewPos ','velocityScrewNeg ', & 'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]),pInt) * prm%totalNslip !< other dependent state variables that are not updated by microstructure sizeDeltaState = sizeDotState call material_allocatePlasticState(p,NofMyPhase,sizeState,sizeDotState,sizeDeltaState, & prm%totalNslip,0_pInt,0_pInt) plasticState(p)%nonlocal = .true. plasticState(p)%offsetDeltaState = 0_pInt ! ToDo: state structure does not follow convention plasticState(p)%sizePostResults = sum(plastic_nonlocal_sizePostResult(:,phase_plasticityInstance(p))) totalNslip(phase_plasticityInstance(p)) = prm%totalNslip ! ToDo: Not really sure if this large number of mostly overlapping pointers is useful stt%rho => plasticState(p)%state (0_pInt*prm%totalNslip+1_pInt:10_pInt*prm%totalNslip,:) dot%rho => plasticState(p)%dotState (0_pInt*prm%totalNslip+1_pInt:10_pInt*prm%totalNslip,:) del%rho => plasticState(p)%deltaState (0_pInt*prm%totalNslip+1_pInt:10_pInt*prm%totalNslip,:) plasticState(p)%aTolState(1:10_pInt*prm%totalNslip) = prm%aTolRho stt%rhoSglEdge => plasticState(p)%state (0_pInt*prm%totalNslip+1_pInt:06_pInt*prm%totalNslip:2*prm%totalNslip,:) stt%rhoSglScrew => plasticState(p)%state (2_pInt*prm%totalNslip+1_pInt:08_pInt*prm%totalNslip:2*prm%totalNslip,:) stt%rhoSgl => plasticState(p)%state (0_pInt*prm%totalNslip+1_pInt: 8_pInt*prm%totalNslip,:) dot%rhoSgl => plasticState(p)%dotState (0_pInt*prm%totalNslip+1_pInt: 8_pInt*prm%totalNslip,:) del%rhoSgl => plasticState(p)%deltaState (0_pInt*prm%totalNslip+1_pInt: 8_pInt*prm%totalNslip,:) stt%rhoSglMobile => plasticState(p)%state (0_pInt*prm%totalNslip+1_pInt: 4_pInt*prm%totalNslip,:) dot%rhoSglMobile => plasticState(p)%dotState (0_pInt*prm%totalNslip+1_pInt: 4_pInt*prm%totalNslip,:) del%rhoSglMobile => plasticState(p)%deltaState (0_pInt*prm%totalNslip+1_pInt: 4_pInt*prm%totalNslip,:) stt%rhoSglEdgeMobile => plasticState(p)%state (0_pInt*prm%totalNslip+1_pInt: 2_pInt*prm%totalNslip,:) dot%rhoSglEdgeMobile => plasticState(p)%dotState (0_pInt*prm%totalNslip+1_pInt: 2_pInt*prm%totalNslip,:) del%rhoSglEdgeMobile => plasticState(p)%deltaState (0_pInt*prm%totalNslip+1_pInt: 2_pInt*prm%totalNslip,:) stt%rho_sgl_mob_edg_pos => plasticState(p)%state (0_pInt*prm%totalNslip+1_pInt: 1_pInt*prm%totalNslip,:) dot%rho_sgl_mob_edg_pos => plasticState(p)%dotState (0_pInt*prm%totalNslip+1_pInt: 1_pInt*prm%totalNslip,:) del%rho_sgl_mob_edg_pos => plasticState(p)%deltaState (0_pInt*prm%totalNslip+1_pInt: 1_pInt*prm%totalNslip,:) stt%rho_sgl_mob_edg_neg => plasticState(p)%state (1_pInt*prm%totalNslip+1_pInt: 2_pInt*prm%totalNslip,:) dot%rho_sgl_mob_edg_neg => plasticState(p)%dotState (1_pInt*prm%totalNslip+1_pInt: 2_pInt*prm%totalNslip,:) del%rho_sgl_mob_edg_neg => plasticState(p)%deltaState (1_pInt*prm%totalNslip+1_pInt: 2_pInt*prm%totalNslip,:) stt%rhoSglScrewMobile => plasticState(p)%state (2_pInt*prm%totalNslip+1_pInt: 4_pInt*prm%totalNslip,:) dot%rhoSglScrewMobile => plasticState(p)%dotState (2_pInt*prm%totalNslip+1_pInt: 4_pInt*prm%totalNslip,:) del%rhoSglScrewMobile => plasticState(p)%deltaState (2_pInt*prm%totalNslip+1_pInt: 4_pInt*prm%totalNslip,:) stt%rho_sgl_mob_scr_pos => plasticState(p)%state (2_pInt*prm%totalNslip+1_pInt: 3_pInt*prm%totalNslip,:) dot%rho_sgl_mob_scr_pos => plasticState(p)%dotState (2_pInt*prm%totalNslip+1_pInt: 3_pInt*prm%totalNslip,:) del%rho_sgl_mob_scr_pos => plasticState(p)%deltaState (2_pInt*prm%totalNslip+1_pInt: 3_pInt*prm%totalNslip,:) stt%rho_sgl_mob_scr_neg => plasticState(p)%state (3_pInt*prm%totalNslip+1_pInt: 4_pInt*prm%totalNslip,:) dot%rho_sgl_mob_scr_neg => plasticState(p)%dotState (3_pInt*prm%totalNslip+1_pInt: 4_pInt*prm%totalNslip,:) del%rho_sgl_mob_scr_neg => plasticState(p)%deltaState (3_pInt*prm%totalNslip+1_pInt: 4_pInt*prm%totalNslip,:) stt%rhoSglImmobile => plasticState(p)%state (4_pInt*prm%totalNslip+1_pInt: 8_pInt*prm%totalNslip,:) dot%rhoSglImmobile => plasticState(p)%dotState (4_pInt*prm%totalNslip+1_pInt: 8_pInt*prm%totalNslip,:) del%rhoSglImmobile => plasticState(p)%deltaState (4_pInt*prm%totalNslip+1_pInt: 8_pInt*prm%totalNslip,:) stt%rhoSglEdgeImmobile => plasticState(p)%state (4_pInt*prm%totalNslip+1_pInt: 6_pInt*prm%totalNslip,:) dot%rhoSglEdgeImmobile => plasticState(p)%dotState (4_pInt*prm%totalNslip+1_pInt: 6_pInt*prm%totalNslip,:) del%rhoSglEdgeImmobile => plasticState(p)%deltaState (4_pInt*prm%totalNslip+1_pInt: 6_pInt*prm%totalNslip,:) stt%rho_sgl_imm_edg_pos => plasticState(p)%state (4_pInt*prm%totalNslip+1_pInt: 5_pInt*prm%totalNslip,:) dot%rho_sgl_imm_edg_pos => plasticState(p)%dotState (4_pInt*prm%totalNslip+1_pInt: 5_pInt*prm%totalNslip,:) del%rho_sgl_imm_edg_pos => plasticState(p)%deltaState (4_pInt*prm%totalNslip+1_pInt: 5_pInt*prm%totalNslip,:) stt%rho_sgl_imm_edg_neg => plasticState(p)%state (5_pInt*prm%totalNslip+1_pInt: 6_pInt*prm%totalNslip,:) dot%rho_sgl_imm_edg_neg => plasticState(p)%dotState (5_pInt*prm%totalNslip+1_pInt: 6_pInt*prm%totalNslip,:) del%rho_sgl_imm_edg_neg => plasticState(p)%deltaState (5_pInt*prm%totalNslip+1_pInt: 6_pInt*prm%totalNslip,:) stt%rhoSglScrewImmobile => plasticState(p)%state (6_pInt*prm%totalNslip+1_pInt: 8_pInt*prm%totalNslip,:) dot%rhoSglScrewImmobile => plasticState(p)%dotState (6_pInt*prm%totalNslip+1_pInt: 8_pInt*prm%totalNslip,:) del%rhoSglScrewImmobile => plasticState(p)%deltaState (6_pInt*prm%totalNslip+1_pInt: 8_pInt*prm%totalNslip,:) stt%rho_sgl_imm_scr_pos => plasticState(p)%state (6_pInt*prm%totalNslip+1_pInt: 7_pInt*prm%totalNslip,:) dot%rho_sgl_imm_scr_pos => plasticState(p)%dotState(6_pInt*prm%totalNslip+1_pInt: 7_pInt*prm%totalNslip,:) del%rho_sgl_imm_scr_pos => plasticState(p)%deltaState(6_pInt*prm%totalNslip+1_pInt: 7_pInt*prm%totalNslip,:) stt%rho_sgl_imm_scr_neg => plasticState(p)%state (7_pInt*prm%totalNslip+1_pInt: 8_pInt*prm%totalNslip,:) dot%rho_sgl_imm_scr_neg => plasticState(p)%dotState(7_pInt*prm%totalNslip+1_pInt: 8_pInt*prm%totalNslip,:) del%rho_sgl_imm_scr_neg => plasticState(p)%deltaState(7_pInt*prm%totalNslip+1_pInt: 8_pInt*prm%totalNslip,:) stt%rhoDip => plasticState(p)%state (8_pInt*prm%totalNslip+1_pInt:10_pInt*prm%totalNslip,:) dot%rhoDip => plasticState(p)%dotState (8_pInt*prm%totalNslip+1_pInt:10_pInt*prm%totalNslip,:) del%rhoDip => plasticState(p)%deltaState (8_pInt*prm%totalNslip+1_pInt:10_pInt*prm%totalNslip,:) stt%rho_dip_edg => plasticState(p)%state (8_pInt*prm%totalNslip+1_pInt: 9_pInt*prm%totalNslip,:) dot%rho_dip_edg => plasticState(p)%dotState (8_pInt*prm%totalNslip+1_pInt: 9_pInt*prm%totalNslip,:) del%rho_dip_edg => plasticState(p)%deltaState (8_pInt*prm%totalNslip+1_pInt: 9_pInt*prm%totalNslip,:) stt%rho_dip_scr => plasticState(p)%state (9_pInt*prm%totalNslip+1_pInt:10_pInt*prm%totalNslip,:) dot%rho_dip_scr => plasticState(p)%dotState (9_pInt*prm%totalNslip+1_pInt:10_pInt*prm%totalNslip,:) del%rho_dip_scr => plasticState(p)%deltaState (9_pInt*prm%totalNslip+1_pInt:10_pInt*prm%totalNslip,:) stt%accumulatedshear => plasticState(p)%state (10_pInt*prm%totalNslip + 1_pInt:11_pInt*prm%totalNslip ,1:NofMyPhase) dot%accumulatedshear => plasticState(p)%dotState (10_pInt*prm%totalNslip + 1_pInt:11_pInt*prm%totalNslip ,1:NofMyPhase) del%accumulatedshear => plasticState(p)%deltaState (10_pInt*prm%totalNslip + 1_pInt:11_pInt*prm%totalNslip ,1:NofMyPhase) plasticState(p)%aTolState(10_pInt*prm%totalNslip + 1_pInt:11_pInt*prm%totalNslip ) = prm%aTolShear plasticState(p)%slipRate => plasticState(p)%dotState(10_pInt*prm%totalNslip + 1_pInt:11_pInt*prm%totalNslip ,1:NofMyPhase) plasticState(p)%accumulatedSlip => plasticState(p)%state (10_pInt*prm%totalNslip + 1_pInt:11_pInt*prm%totalNslip ,1:NofMyPhase) stt%rho_forest => plasticState(p)%state (11_pInt*prm%totalNslip + 1_pInt:12_pInt*prm%totalNslip ,1:NofMyPhase) allocate(dst%tau_Threshold(prm%totalNslip,NofMyPhase),source=0.0_pReal) allocate(dst%tau_Back(prm%totalNslip,NofMyPhase),source=0.0_pReal) allocate(res%rhoDotFlux(prm%totalNslip,8,NofMyPhase),source=0.0_pReal) allocate(res%rhoDotMultiplication(prm%totalNslip,2,NofMyPhase),source=0.0_pReal) allocate(res%rhoDotSingle2DipoleGlide(prm%totalNslip,2,NofMyPhase),source=0.0_pReal) allocate(res%rhoDotAthermalAnnihilation(prm%totalNslip,2,NofMyPhase),source=0.0_pReal) allocate(res%rhoDotThermalAnnihilation(prm%totalNslip,2,NofMyPhase),source=0.0_pReal) allocate(res%rhoDotEdgeJogs(prm%totalNslip,NofMyPhase),source=0.0_pReal) end associate if (NofMyPhase > 0_pInt) call stateInit(p,NofMyPhase) plasticState(p)%state0 = plasticState(p)%state enddo ! BEGIN DEPRECATED---------------------------------------------------------------------------------- allocate(iRhoU(maxval(totalNslip),4,maxNinstances), source=0_pInt) allocate(iRhoB(maxval(totalNslip),4,maxNinstances), source=0_pInt) allocate(iRhoD(maxval(totalNslip),2,maxNinstances), source=0_pInt) allocate(iV(maxval(totalNslip),4,maxNinstances), source=0_pInt) allocate(iD(maxval(totalNslip),2,maxNinstances), source=0_pInt) allocate(iRhoF(maxval(totalNslip),maxNinstances), source=0_pInt) ! END DEPRECATED------------------------------------------------------------------------------------ allocate(compatibility(2,maxval(totalNslip),maxval(totalNslip),theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems), & source=0.0_pReal) initializeInstances: do p = 1_pInt, size(phase_plasticity) NofMyPhase=count(material_phase==p) myPhase2: if (phase_plasticity(p) == PLASTICITY_NONLOCAL_ID) then !*** determine indices to state array l = 0_pInt do t = 1_pInt,4_pInt do s = 1_pInt,param(phase_plasticityInstance(p))%totalNslip l = l + 1_pInt iRhoU(s,t,phase_plasticityInstance(p)) = l enddo enddo do t = 1_pInt,4_pInt do s = 1_pInt,param(phase_plasticityInstance(p))%totalNslip l = l + 1_pInt iRhoB(s,t,phase_plasticityInstance(p)) = l enddo enddo do c = 1_pInt,2_pInt do s = 1_pInt,param(phase_plasticityInstance(p))%totalNslip l = l + 1_pInt iRhoD(s,c,phase_plasticityInstance(p)) = l enddo enddo l = l + param(phase_plasticityInstance(p))%totalNslip ! shear(rates) do s = 1_pInt,param(phase_plasticityInstance(p))%totalNslip l = l + 1_pInt iRhoF(s,phase_plasticityInstance(p)) = l enddo do t = 1_pInt,4_pInt do s = 1_pInt,param(phase_plasticityInstance(p))%totalNslip l = l + 1_pInt iV(s,t,phase_plasticityInstance(p)) = l enddo enddo do c = 1_pInt,2_pInt do s = 1_pInt,param(phase_plasticityInstance(p))%totalNslip l = l + 1_pInt iD(s,c,phase_plasticityInstance(p)) = l enddo enddo if (iD(param(phase_plasticityInstance(p))%totalNslip,2,phase_plasticityInstance(p)) /= plasticState(p)%sizeState) & ! check if last index is equal to size of state call IO_error(0_pInt, ext_msg = 'state indices not properly set ('//PLASTICITY_NONLOCAL_label//')') endif myPhase2 enddo initializeInstances contains subroutine stateInit(phase,NofMyPhase) use math, only: & math_sampleGaussVar use mesh, only: & theMesh, & mesh_ipVolume use material, only: & material_phase, & phase_plasticityInstance, & phasememberAt implicit none integer(pInt),intent(in) ::& phase, & NofMyPhase integer(pInt) :: & e, & i, & f, & from, & upto, & s, & instance, & phasemember real(pReal), dimension(2) :: & noise, & rnd real(pReal) :: & meanDensity, & totalVolume, & densityBinning, & 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 if (prm%rhoSglRandom > 0.0_pReal) then ! get the total volume of the instance do e = 1_pInt,theMesh%nElems do i = 1_pInt,theMesh%elem%nIPs if (material_phase(1,i,e) == phase) volume(phasememberAt(1,i,e)) = mesh_ipVolume(i,e) enddo enddo 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) call random_number(rnd) phasemember = nint(rnd(1)*real(NofMyPhase,pReal) + 0.5_pReal,pInt) s = nint(rnd(2)*real(prm%totalNslip,pReal)*4.0_pReal + 0.5_pReal,pInt) meanDensity = meanDensity + densityBinning * volume(phasemember) / totalVolume stt%rhoSglMobile(s,phasemember) = densityBinning enddo ! homogeneous distribution of density with some noise else do e = 1_pInt, NofMyPhase do f = 1_pInt,size(prm%Nslip,1) from = 1_pInt + sum(prm%Nslip(1:f-1_pInt)) upto = sum(prm%Nslip(1:f)) do s = from,upto noise = [math_sampleGaussVar(0.0_pReal, prm%rhoSglScatter), & math_sampleGaussVar(0.0_pReal, prm%rhoSglScatter)] stt%rho_sgl_mob_edg_pos(s,e) = prm%rhoSglEdgePos0(f) + noise(1) stt%rho_sgl_mob_edg_neg(s,e) = prm%rhoSglEdgeNeg0(f) + noise(1) stt%rho_sgl_mob_scr_pos(s,e) = prm%rhoSglScrewPos0(f) + noise(2) stt%rho_sgl_mob_scr_neg(s,e) = prm%rhoSglScrewNeg0(f) + noise(2) enddo stt%rho_dip_edg(from:upto,e) = prm%rhoDipEdge0(f) stt%rho_dip_scr(from:upto,e) = prm%rhoDipScrew0(f) enddo enddo endif end associate end subroutine stateInit end subroutine plastic_nonlocal_init !-------------------------------------------------------------------------------------------------- !> @brief calculates quantities characterizing the microstructure !-------------------------------------------------------------------------------------------------- subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) use prec, only: & dEq0 use IO, only: & IO_error use math, only: & pi, & math_mul33x3, & math_mul3x3, & math_inv33 #ifdef DEBUG use debug, only: & debug_level, & debug_constitutive, & debug_levelExtensive, & debug_levelSelective, & debug_i, & debug_e #endif use mesh, only: & theMesh, & mesh_ipNeighborhood, & mesh_ipCoordinates, & mesh_ipVolume, & mesh_ipAreaNormal, & mesh_ipArea use material, only: & material_phase, & phase_localPlasticity, & phaseAt, phasememberAt, & phase_plasticityInstance use lattice, only: & LATTICE_bcc_ID, & LATTICE_fcc_ID, & lattice_structure implicit none integer(pInt), intent(in) :: ip, & ! current integration point el ! current element real(pReal), dimension(3,3), intent(in) :: & Fe, & ! elastic deformation gradient Fp ! elastic deformation gradient integer(pInt) :: & ph, & !< phase of, & !< offset no !< nieghbor offset integer(pInt) ns, neighbor_el, & ! element number of neighboring material point neighbor_ip, & ! integration point of neighboring material point instance, & ! my instance of this plasticity neighbor_instance, & ! instance of this plasticity of neighboring material point c, & ! index of dilsocation character (edge, screw) s, & ! slip system index t, & ! index of dilsocation type (e+, e-, s+, s-, used e+, used e-, used s+, used s-) dir, & n, & nRealNeighbors ! number of really existing neighbors integer(pInt), dimension(2) :: neighbors real(pReal) FVsize, & correction real(pReal), dimension(2) :: rhoExcessGradient, & rhoExcessGradient_over_rho, & rhoTotal real(pReal), dimension(3) :: rhoExcessDifferences, & normal_latticeConf real(pReal), dimension(3,3) :: invFe, & ! inverse of elastic deformation gradient invFp, & ! inverse of plastic deformation gradient connections, & invConnections real(pReal), dimension(3,theMesh%elem%nIPneighbors) :: & connection_latticeConf real(pReal), dimension(2,totalNslip(phase_plasticityInstance(material_phase(1_pInt,ip,el)))) :: & rhoExcess real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el)))) :: & rho_edg_delta, & rho_scr_delta real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),10) :: & rho, & rho_neighbor real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1_pInt,ip,el))),8) :: & rhoSgl ! single dislocation density (edge+, edge-, screw+, screw-, used edge+, used edge-, used screw+, used screw-) real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1_pInt,ip,el))), & totalNslip(phase_plasticityInstance(material_phase(1_pInt,ip,el)))) :: & myInteractionMatrix ! corrected slip interaction matrix real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),theMesh%elem%nIPneighbors) :: & rho_edg_delta_neighbor, & rho_scr_delta_neighbor, & rho_edg_sum_neighbor, & rho_scr_sum_neighbor real(pReal), dimension(2,maxval(totalNslip),theMesh%elem%nIPneighbors) :: & neighbor_rhoExcess, & ! excess density at neighboring material point neighbor_rhoTotal ! total density at neighboring material point real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(1_pInt,ip,el))),2) :: & m ! direction of dislocation motion ph = phaseAt(1,ip,el) of = 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) rhoSgl = rho(:,sgl) stt%rho_forest(:,of) = matmul(prm%forestProjection_Edge, sum(abs(rho(:,edg)),2)) & + matmul(prm%forestProjection_Screw,sum(abs(rho(:,scr)),2)) !*** calculate the threshold shear stress for dislocation slip !*** coefficients are corrected for the line tension effect !*** (see Kubin,Devincre,Hoc; 2008; Modeling dislocation storage rates and mean free paths in face-centered cubic crystals) if (lattice_structure(ph) == LATTICE_bcc_ID .or. lattice_structure(ph) == LATTICE_fcc_ID) then ! only fcc and bcc do s = 1_pInt,ns correction = ( 1.0_pReal - prm%linetensionEffect & + 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) enddo else myInteractionMatrix = prm%interactionSlipSlip endif forall (s = 1_pInt:ns) & dst%tau_threshold(s,of) = prm%mu * prm%burgers(s) & * sqrt(dot_product(sum(abs(rho),2), myInteractionMatrix(1:ns,s))) !*** calculate the dislocation stress of the neighboring excess dislocation densities !*** zero for material points of local plasticity dst%tau_back(:,of) = 0.0_pReal !################################################################################################# !################################################################################################# ! ToDo: MD: this is most likely only correct for F_i = I !################################################################################################# !################################################################################################# 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 = mesh_ipVolume(ip,el) ** (1.0_pReal/3.0_pReal) !* loop through my neighborhood and get the connection vectors (in lattice frame) and the excess densities nRealNeighbors = 0 neighbor_rhoTotal = 0.0_pReal do n = 1_pInt,theMesh%elem%nIPneighbors neighbor_el = mesh_ipNeighborhood(1,n,ip,el) neighbor_ip = mesh_ipNeighborhood(2,n,ip,el) no = phasememberAt(1,neighbor_ip,neighbor_el) if (neighbor_el > 0 .and. neighbor_ip > 0) then neighbor_instance = phase_plasticityInstance(material_phase(1,neighbor_ip,neighbor_el)) if (neighbor_instance == instance) then nRealNeighbors = nRealNeighbors + 1 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) = & math_mul33x3(invFe, mesh_ipCoordinates(1:3,neighbor_ip,neighbor_el) & - mesh_ipCoordinates(1:3,ip,el)) normal_latticeConf = math_mul33x3(transpose(invFp), mesh_ipAreaNormal(1:3,n,ip,el)) if (math_mul3x3(normal_latticeConf,connection_latticeConf(1:3,n)) < 0.0_pReal) & ! neighboring connection points in opposite direction to face normal: must be periodic image connection_latticeConf(1:3,n) = normal_latticeConf * mesh_ipVolume(ip,el) & / mesh_ipArea(n,ip,el) ! instead take the surface normal scaled with the diameter of the cell else ! local neighbor or different lattice structure or different constitution instance -> use central values instead connection_latticeConf(1:3,n) = 0.0_pReal rho_edg_delta_neighbor(:,n) = rho_edg_delta rho_scr_delta_neighbor(:,n) = rho_scr_delta endif else ! free surface -> use central values instead connection_latticeConf(1:3,n) = 0.0_pReal rho_edg_delta_neighbor(:,n) = rho_edg_delta 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_pInt,ns ! gradient from interpolation of neighboring excess density ... do c = 1_pInt,2_pInt do dir = 1_pInt,3_pInt neighbors(1) = 2_pInt * dir - 1_pInt neighbors(2) = 2_pInt * dir connections(dir,1:3) = connection_latticeConf(1:3,neighbors(1)) & - connection_latticeConf(1:3,neighbors(2)) rhoExcessDifferences(dir) = neighbor_rhoExcess(c,s,neighbors(1)) & - neighbor_rhoExcess(c,s,neighbors(2)) enddo invConnections = math_inv33(connections) if (all(dEq0(invConnections))) call IO_error(-1_pInt,ext_msg='back stress calculation: inversion error') rhoExcessGradient(c) = math_mul3x3(m(1:3,s,c), math_mul33x3(invConnections,rhoExcessDifferences)) enddo ! ... plus gradient from deads ... do t = 1_pInt,4_pInt c = (t - 1_pInt) / 2_pInt + 1_pInt rhoExcessGradient(c) = rhoExcessGradient(c) + rhoSgl(s,t+4_pInt) / FVsize enddo ! ... normalized with the total density ... rhoTotal(1) = (sum(abs(rho(s,edg))) + sum(neighbor_rhoTotal(1,s,:))) / real(1_pInt + nRealNeighbors,pReal) rhoTotal(2) = (sum(abs(rho(s,scr))) + sum(neighbor_rhoTotal(2,s,:))) / real(1_pInt + nRealNeighbors,pReal) 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 #ifdef DEBUG if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt & .and. ((debug_e == el .and. debug_i == ip)& .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then write(6,'(/,a,i8,1x,i2,1x,i1,/)') '<< CONST >> nonlocal_microstructure at el ip ',el,ip write(6,'(a,/,12x,12(e10.3,1x))') '<< CONST >> rhoForest', stt%rho_forest(:,of) write(6,'(a,/,12x,12(f10.5,1x))') '<< CONST >> tauThreshold / MPa', dst%tau_threshold(:,of)*1e-6 write(6,'(a,/,12x,12(f10.5,1x),/)') '<< CONST >> tauBack / MPa', dst%tau_back(:,of)*1e-6 endif #endif end associate end subroutine plastic_nonlocal_dependentState !-------------------------------------------------------------------------------------------------- !> @brief calculates kinetics !-------------------------------------------------------------------------------------------------- subroutine plastic_nonlocal_kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, & tauThreshold, c, Temperature, instance, of) implicit none integer(pInt), intent(in) :: c, & !< dislocation character (1:edge, 2:screw) instance, of real(pReal), intent(in) :: Temperature !< temperature real(pReal), dimension(param(instance)%totalNslip), & intent(in) :: 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(pInt) :: ns, & !< short notation for the total number of active slip systems s !< index of my current slip system real(pReal) tauRel_P, & tauRel_S, & tauEff, & !< effective shear stress tPeierls, & !< waiting time in front of a peierls barriers tSolidSolution, & !< waiting time in front of a solid solution obstacle vViscous, & !< viscous glide velocity dtPeierls_dtau, & !< derivative with respect to resolved shear stress dtSolidSolution_dtau, & !< derivative with respect to resolved shear stress meanfreepath_S, & !< mean free travel distance for dislocations between two solid solution obstacles meanfreepath_P, & !< mean free travel distance for dislocations between two Peierls barriers jumpWidth_P, & !< depth of activated area jumpWidth_S, & !< depth of activated area activationLength_P, & !< length of activated dislocation line activationLength_S, & !< length of activated dislocation line activationVolume_P, & !< volume that needs to be activated to overcome barrier activationVolume_S, & !< volume that needs to be activated to overcome barrier activationEnergy_P, & !< energy that is needed to overcome barrier activationEnergy_S, & !< energy that is needed to overcome barrier 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_pInt,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) activationLength_P = prm%doublekinkwidth *prm%burgers(s) activationVolume_P = activationLength_P * jumpWidth_P * prm%burgers(s) criticalStress_P = prm%peierlsStress(s,c) activationEnergy_P = criticalStress_P * activationVolume_P tauRel_P = min(1.0_pReal, tauEff / criticalStress_P) ! ensure that the activation probability cannot become greater than one tPeierls = 1.0_pReal / prm%fattack & * exp(activationEnergy_P / (KB * Temperature) & * (1.0_pReal - tauRel_P**prm%p)**prm%q) 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) 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) activationLength_S = prm%burgers(s) / sqrt(prm%solidSolutionConcentration) activationVolume_S = activationLength_S * jumpWidth_S * prm%burgers(s) activationEnergy_S = prm%solidSolutionEnergy criticalStress_S = activationEnergy_S / activationVolume_S tauRel_S = min(1.0_pReal, tauEff / criticalStress_S) ! ensure that the activation probability cannot become greater than one tSolidSolution = 1.0_pReal / prm%fattack & * exp(activationEnergy_S / (KB * Temperature) & * (1.0_pReal - tauRel_S**prm%p)**prm%q) if (tauEff < criticalStress_S) then 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) 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. !* 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 endif enddo endif #ifdef DEBUGTODO write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> tauThreshold / MPa', tauThreshold * 1e-6_pReal write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> tau / MPa', tau * 1e-6_pReal write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> tauNS / MPa', tauNS * 1e-6_pReal write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> v / mm/s', v * 1e3 write(6,'(a,/,12x,12(e12.5,1x))') '<< CONST >> dv_dtau', dv_dtau write(6,'(a,/,12x,12(e12.5,1x))') '<< CONST >> dv_dtauNS', dv_dtauNS endif #endif end associate end subroutine plastic_nonlocal_kinetics !-------------------------------------------------------------------------------------------------- !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dMp, & Mp, Temperature, volume, ip, el) use math, only: & math_mul33xx33 use material, only: & material_phase, & plasticState, & phaseAt, phasememberAt,& phase_plasticityInstance implicit none integer(pInt), intent(in) :: ip, & !< current integration point el !< current element number real(pReal), intent(in) :: Temperature, & !< temperature volume !< volume of the materialpoint real(pReal), dimension(3,3), intent(in) :: Mp 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 Tstar (9x9 matrix) integer(pInt) instance, & !< current instance of this plasticity ns, & !< short notation for the total number of active slip systems i, & j, & k, & l, & ph, & !phase number of, & !offset t, & !< dislocation type s !< index of my current slip system real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1_pInt,ip,el))),8) :: & rhoSgl !< single dislocation densities (including blocked) real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1_pInt,ip,el))),10) :: & rho real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1_pInt,ip,el))),4) :: & v, & !< velocity tauNS, & !< resolved shear stress including non Schmid and backstress terms dv_dtau, & !< velocity derivative with respect to the shear stress dv_dtauNS !< velocity derivative with respect to the shear stress real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1_pInt,ip,el)))) :: & tau, & !< resolved shear stress including backstress terms gdotTotal !< shear rate !*** shortcut for mapping ph = phaseAt(1_pInt,ip,el) of = phasememberAt(1_pInt,ip,el) instance = phase_plasticityInstance(ph) associate(prm => param(instance),dst=>microstructure(instance)) ns = prm%totalNslip !*** 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_pInt,ns tau(s) = math_mul33xx33(Mp, prm%Schmid(1:3,1:3,s)) tauNS(s,1) = tau(s) tauNS(s,2) = tau(s) if (tau(s) > 0.0_pReal) then tauNS(s,3) = math_mul33xx33(Mp, +prm%nonSchmid_pos(1:3,1:3,s)) tauNS(s,4) = math_mul33xx33(Mp, -prm%nonSchmid_neg(1:3,1:3,s)) else tauNS(s,3) = math_mul33xx33(Mp, +prm%nonSchmid_neg(1:3,1:3,s)) tauNS(s,4) = math_mul33xx33(Mp, -prm%nonSchmid_pos(1:3,1:3,s)) endif enddo forall (t = 1_pInt:4_pInt) & 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 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_Threshold(1:ns,of), & 1_pInt, 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_pInt) then ! no non-Schmid contributions forall(t = 3_pInt:4_pInt) v(1:ns,t) = v(1:ns,1) dv_dtau(1:ns,t) = dv_dtau(1:ns,1) dv_dtauNS(1:ns,t) = dv_dtauNS(1:ns,1) endforall else ! take non-Schmid contributions into account do t = 3_pInt,4_pInt call plastic_nonlocal_kinetics(v(1:ns,t), dv_dtau(1:ns,t), dv_dtauNS(1:ns,t), & tau(1:ns), tauNS(1:ns,t), dst%tau_Threshold(1:ns,of), & 2_pInt , Temperature, instance, of) enddo endif !*** store velocity in state forall (t = 1_pInt:4_pInt) & plasticState(ph)%state(iV(1:ns,t,instance),of) = v(1:ns,t) !*** Bauschinger effect forall (s = 1_pInt:ns, t = 5_pInt:8_pInt, rhoSgl(s,t) * v(s,t-4_pInt) < 0.0_pReal) & rhoSgl(s,t-4_pInt) = rhoSgl(s,t-4_pInt) + abs(rhoSgl(s,t)) !*** Calculation of Lp and its tangent 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_pInt,ns Lp = Lp + gdotTotal(s) * prm%Schmid(1:3,1:3,s) forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt,k=1_pInt:3_pInt,l=1_pInt:3_pInt) & dLp_dMp(i,j,k,l) = dLp_dMp(i,j,k,l) & + 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_neg(k,l,s) * rhoSgl(s,4) * dv_dtauNS(s,4)) * prm%burgers(s) enddo end associate end subroutine plastic_nonlocal_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief (instantaneous) incremental change of microstructure !-------------------------------------------------------------------------------------------------- subroutine plastic_nonlocal_deltaState(Mp,ip,el) use prec, only: & dNeq0 use debug, only: debug_level, & debug_constitutive, & debug_levelBasic, & debug_levelExtensive, & debug_levelSelective, & debug_i, & debug_e use math, only: PI, & math_mul33xx33 use mesh, only: mesh_ipVolume use material, only: material_phase, & plasticState, & phaseAt, phasememberAt, & phase_plasticityInstance implicit none integer(pInt), intent(in) :: ip, & ! current grain number el ! current element number real(pReal), dimension(3,3), intent(in) :: Mp !< MandelStress integer(pInt) :: & ph, & !< phase of !< offset integer(pInt) ::instance, & ! current instance of this plasticity ns, & ! short notation for the total number of active slip systems c, & ! character of dislocation t, & ! type of dislocation s ! index of my current slip system real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),10) :: & deltaRho, & ! density increment deltaRhoRemobilization, & ! density increment by remobilization deltaRhoDipole2SingleStress ! density increment by dipole dissociation (by stress change) real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),10) :: & rho ! current dislocation densities real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),8) :: & rhoSgl ! current single dislocation densities (positive/negative screw and edge without dipoles) real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),4) :: & v ! dislocation glide velocity real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el)))) :: & tau ! current resolved shear stress real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),2) :: & rhoDip, & ! current dipole dislocation densities (screw and edge dipoles) dLower, & ! minimum stable dipole distance for edges and screws dUpper, & ! current maximum stable dipole distance for edges and screws dUpperOld, & ! old maximum stable dipole distance for edges and screws deltaDUpper ! change in maximum stable dipole distance for edges and screws #ifdef DEBUG if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt & .and. ((debug_e == el .and. debug_i == ip)& .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) & write(6,'(/,a,i8,1x,i2,1x,i1,/)') '<< CONST >> nonlocal_deltaState at el ip ',el,ip #endif ph = phaseAt(1,ip,el) of = phasememberAt(1,ip,el) instance = phase_plasticityInstance(ph) associate(prm => param(instance),dst => microstructure(instance)) ns = totalNslip(instance) !*** shortcut to state variables forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) v(s,t) = plasticState(ph)%state(iV(s,t,instance),of) endforall forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,instance),of) endforall rho = getRho(instance,of,ip,el) rhoSgl = rho(:,sgl) rhoDip = rho(:,dip) !**************************************************************************** !*** dislocation remobilization (bauschinger effect) deltaRhoRemobilization = 0.0_pReal do t = 1_pInt,4_pInt do s = 1_pInt,ns if (rhoSgl(s,t+4_pInt) * v(s,t) < 0.0_pReal) then deltaRhoRemobilization(s,t) = abs(rhoSgl(s,t+4_pInt)) rhoSgl(s,t) = rhoSgl(s,t) + abs(rhoSgl(s,t+4_pInt)) deltaRhoRemobilization(s,t+4_pInt) = - rhoSgl(s,t+4_pInt) rhoSgl(s,t+4_pInt) = 0.0_pReal endif enddo enddo !**************************************************************************** !*** calculate dipole formation and dissociation by stress change !*** calculate limits for stable dipole height do s = 1_pInt,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 dLower = prm%minDipoleHeight(1:ns,1:2) 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)) do c = 1, 2 where(dNeq0(sqrt(rhoSgl(1:ns,2*c-1)+rhoSgl(1:ns,2*c)+abs(rhoSgl(1:ns,2*c+3))& +abs(rhoSgl(1:ns,2*c+4))+rhoDip(1:ns,c)))) & dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) & + abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), & dUpper(1:ns,c)) enddo dUpper = max(dUpper,dLower) deltaDUpper = dUpper - dUpperOld !*** dissociation by stress increase deltaRhoDipole2SingleStress = 0.0_pReal forall (c=1_pInt:2_pInt, s=1_pInt:ns, deltaDUpper(s,c) < 0.0_pReal .and. & dNeq0(dUpperOld(s,c) - dLower(s,c))) & deltaRhoDipole2SingleStress(s,8_pInt+c) = rhoDip(s,c) * deltaDUpper(s,c) & / (dUpperOld(s,c) - dLower(s,c)) forall (t=1_pInt:4_pInt) & deltaRhoDipole2SingleStress(1_pInt:ns,t) = -0.5_pReal & * deltaRhoDipole2SingleStress(1_pInt:ns,(t-1_pInt)/2_pInt+9_pInt) !*** store new maximum dipole height in state forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) & plasticState(ph)%state(iD(s,c,instance),of) = dUpper(s,c) !**************************************************************************** !*** assign the changes in the dislocation densities to deltaState deltaRho = deltaRhoRemobilization & + deltaRhoDipole2SingleStress plasticState(ph)%deltaState(:,of) = 0.0_pReal forall (s = 1:ns, t = 1_pInt:4_pInt) plasticState(ph)%deltaState(iRhoU(s,t,instance),of)= deltaRho(s,t) plasticState(ph)%deltaState(iRhoB(s,t,instance),of) = deltaRho(s,t+4_pInt) endforall forall (s = 1:ns, c = 1_pInt:2_pInt) & plasticState(ph)%deltaState(iRhoD(s,c,instance),of) = deltaRho(s,c+8_pInt) #ifdef DEBUG if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt & .and. ((debug_e == el .and. debug_i == ip)& .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then write(6,'(a,/,8(12x,12(e12.5,1x),/))') '<< CONST >> dislocation remobilization', deltaRhoRemobilization(1:ns,1:8) write(6,'(a,/,10(12x,12(e12.5,1x),/),/)') '<< CONST >> dipole dissociation by stress increase', deltaRhoDipole2SingleStress endif #endif end associate end subroutine plastic_nonlocal_deltaState !--------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !--------------------------------------------------------------------------------------------------- subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & timestep,ip,el) use, intrinsic :: & IEEE_arithmetic use prec, only: dNeq0, & dNeq, & dEq0 use IO, only: IO_error #ifdef DEBUG use debug, only: debug_level, & debug_constitutive, & debug_levelBasic, & debug_levelExtensive, & debug_levelSelective, & debug_i, & debug_e #endif use math, only: math_mul3x3, & #ifdef __PGI norm2, & #endif math_mul33x3, & math_mul33xx33, & math_mul33x33, & math_inv33, & math_det33, & pi use mesh, only: theMesh, & mesh_ipNeighborhood, & mesh_ipVolume, & mesh_ipArea, & mesh_ipAreaNormal use material, only: homogenization_maxNgrains, & material_phase, & phase_plasticityInstance, & phase_localPlasticity, & plasticState, & phaseAt, phasememberAt, & phase_plasticity ,& PLASTICITY_NONLOCAL_ID use lattice, only: lattice_structure, & LATTICE_bcc_ID, & LATTICE_fcc_ID implicit none !*** input variables integer(pInt), intent(in) :: ip, & !< current integration point el !< current element number real(pReal), intent(in) :: Temperature, & !< temperature timestep !< substepped crystallite time increment real(pReal), dimension(3,3), intent(in) :: Mp !< MandelStress real(pReal), dimension(3,3,homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%nElems), intent(in) :: & Fe, & !< elastic deformation gradient Fp !< plastic deformation gradient !*** local variables integer(pInt) :: ph, & instance, & !< current instance of this plasticity neighbor_instance, & !< instance of my neighbor's plasticity ns, & !< short notation for the total number of active slip systems c, & !< character of dislocation n, & !< index of my current neighbor neighbor_el, & !< element number of my neighbor neighbor_ip, & !< integration point of my neighbor neighbor_n, & !< neighbor index pointing to me when looking from my neighbor opposite_neighbor, & !< index of my opposite neighbor opposite_ip, & !< ip of my opposite neighbor opposite_el, & !< element index of my opposite neighbor opposite_n, & !< neighbor index pointing to me when looking from my opposite neighbor t, & !< type of dislocation o,& !< offset shortcut no,& !< neighbour offset shortcut p,& !< phase shortcut np,& !< neighbour phase shortcut topp, & !< type of dislocation with opposite sign to t s !< index of my current slip system real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1_pInt,ip,el))),10) :: & rho, & rhoDot, & !< density evolution rhoDotMultiplication, & !< density evolution by multiplication rhoDotFlux, & !< density evolution by flux rhoDotSingle2DipoleGlide, & !< density evolution by dipole formation (by glide) rhoDotAthermalAnnihilation, & !< density evolution by athermal annihilation rhoDotThermalAnnihilation !< density evolution by thermal annihilation real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1_pInt,ip,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) my_rhoSgl !< single dislocation densities of central ip (positive/negative screw and edge without dipoles) real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1_pInt,ip,el))),4) :: & v, & !< current dislocation glide velocity my_v, & !< dislocation glide velocity of central ip neighbor_v, & !< dislocation glide velocity of enighboring ip gdot !< shear rates real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1_pInt,ip,el)))) :: & tau, & !< current resolved shear stress vClimb !< climb velocity of edge dipoles real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1_pInt,ip,el))),2) :: & rhoDip, & !< current dipole dislocation densities (screw and edge dipoles) dLower, & !< minimum stable dipole distance for edges and screws dUpper !< current maximum stable dipole distance for edges and screws real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(1_pInt,ip,el))),4) :: & m !< direction of dislocation motion real(pReal), dimension(3,3) :: my_F, & !< my total deformation gradient neighbor_F, & !< total deformation gradient of my neighbor my_Fe, & !< my elastic deformation gradient neighbor_Fe, & !< elastic deformation gradient of my neighbor Favg !< average total deformation gradient of me and my neighbor real(pReal), dimension(3) :: normal_neighbor2me, & !< interface normal pointing from my neighbor to me in neighbor's lattice configuration normal_neighbor2me_defConf, & !< interface normal pointing from my neighbor to me in shared deformed configuration normal_me2neighbor, & !< interface normal pointing from me to my neighbor in my lattice configuration normal_me2neighbor_defConf !< interface normal pointing from me to my neighbor in shared deformed configuration real(pReal) area, & !< area of the current interface transmissivity, & !< overall transmissivity of dislocation flux to neighboring material point lineLength, & !< dislocation line length leaving the current interface selfDiffusion !< self diffusion logical considerEnteringFlux, & considerLeavingFlux p = phaseAt(1,ip,el) o = phasememberAt(1,ip,el) if (timestep <= 0.0_pReal) then ! if illegal timestep... plasticState(p)%dotState = 0.0_pReal ! ...return without doing anything (-> zero dotState) return endif #ifdef DEBUG if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt & .and. ((debug_e == el .and. debug_i == ip)& .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) & write(6,'(/,a,i8,1x,i2,/)') '<< CONST >> nonlocal_dotState at el ip ',el,ip #endif ph = material_phase(1_pInt,ip,el) instance = phase_plasticityInstance(ph) associate(prm => param(instance),dst => microstructure(instance),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_pInt:ns, t = 1_pInt:4_pInt) v(s,t) = plasticState(p)%state(iV (s,t,instance),o) endforall !**************************************************************************** !*** Calculate shear rate forall (t = 1_pInt:4_pInt) & gdot(1_pInt:ns,t) = rhoSgl(1_pInt:ns,t) * prm%burgers(1:ns) * v(1:ns,t) #ifdef DEBUG if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt & .and. ((debug_e == el .and. debug_i == ip)& .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> rho / 1/m^2', rhoSgl, rhoDip write(6,'(a,/,4(12x,12(e12.5,1x),/))') '<< CONST >> gdot / 1/s',gdot endif #endif !**************************************************************************** !*** calculate limits for stable dipole height do s = 1_pInt,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(1:ns,1:2) dUpper(1:ns,1) = prm%mu * prm%burgers(1:ns) & / (8.0_pReal * pi * (1.0_pReal - prm%nu) * abs(tau)) dUpper(1:ns,2) = prm%mu * prm%burgers(1:ns) & / (4.0_pReal * pi * abs(tau)) do c = 1, 2 where(dNeq0(sqrt(rhoSgl(1:ns,2*c-1)+rhoSgl(1:ns,2*c)+abs(rhoSgl(1:ns,2*c+3))& +abs(rhoSgl(1:ns,2*c+4))+rhoDip(1:ns,c)))) & dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) & + abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), & dUpper(1:ns,c)) enddo dUpper = max(dUpper,dLower) !**************************************************************************** !*** calculate dislocation multiplication rhoDotMultiplication = 0.0_pReal if (lattice_structure(ph) == LATTICE_bcc_ID) then ! BCC forall (s = 1:ns, sum(abs(v(s,1:4))) > 0.0_pReal) rhoDotMultiplication(s,1:2) = sum(abs(gdot(s,3:4))) / prm%burgers(s) & ! assuming double-cross-slip of screws to be decisive for multiplication * sqrt(stt%rho_forest(s,o)) / prm%lambda0(s) ! & ! mean free path ! * 2.0_pReal * sum(abs(v(s,3:4))) / sum(abs(v(s,1:4))) ! ratio of screw to overall velocity determines edge generation rhoDotMultiplication(s,3:4) = sum(abs(gdot(s,3:4))) /prm%burgers(s) & ! assuming double-cross-slip of screws to be decisive for multiplication * sqrt(stt%rho_forest(s,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 ! ALL OTHER STRUCTURES 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 !**************************************************************************** !*** calculate dislocation fluxes (only for nonlocal plasticity) rhoDotFlux = 0.0_pReal !? why needed here if (.not. phase_localPlasticity(material_phase(1_pInt,ip,el))) then ! only for nonlocal plasticity !*** 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 & > mesh_ipVolume(ip,el) / maxval(mesh_ipArea(:,ip,el)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here) #ifdef DEBUG if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt) then write(6,'(a,i5,a,i2)') '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip write(6,'(a,e10.3,a,e10.3)') '<< CONST >> velocity is at ', & maxval(abs(v), abs(gdot) > 0.0_pReal & .and. prm%CFLfactor * abs(v) * timestep & > mesh_ipVolume(ip,el) / maxval(mesh_ipArea(:,ip,el))), & ' at a timestep of ',timestep write(6,'(a)') '<< CONST >> enforcing cutback !!!' endif #endif 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_pInt,ip,el) my_F = math_mul33x33(my_Fe, Fp(1:3,1:3,1_pInt,ip,el)) do n = 1_pInt,theMesh%elem%nIPneighbors neighbor_el = mesh_ipNeighborhood(1,n,ip,el) neighbor_ip = mesh_ipNeighborhood(2,n,ip,el) neighbor_n = mesh_ipNeighborhood(3,n,ip,el) np = phaseAt(1,neighbor_ip,neighbor_el) no = phasememberAt(1,neighbor_ip,neighbor_el) opposite_neighbor = n + mod(n,2_pInt) - mod(n+1_pInt,2_pInt) opposite_el = mesh_ipNeighborhood(1,opposite_neighbor,ip,el) opposite_ip = mesh_ipNeighborhood(2,opposite_neighbor,ip,el) opposite_n = mesh_ipNeighborhood(3,opposite_neighbor,ip,el) if (neighbor_n > 0_pInt) then ! if neighbor exists, average deformation gradient neighbor_instance = phase_plasticityInstance(material_phase(1_pInt,neighbor_ip,neighbor_el)) neighbor_Fe = Fe(1:3,1:3,1_pInt,neighbor_ip,neighbor_el) neighbor_F = math_mul33x33(neighbor_Fe, Fp(1:3,1:3,1_pInt,neighbor_ip,neighbor_el)) Favg = 0.5_pReal * (my_F + neighbor_F) else ! if no neighbor, take my value as average Favg = my_F endif !* FLUX FROM MY NEIGHBOR TO ME !* This is only considered, if I have a neighbor of nonlocal plasticity !* (also nonlocal constitutive law with local properties) that is at least a little bit !* compatible. !* If it's not at all compatible, no flux is arriving, because everything is dammed in front of !* my neighbor's interface. !* The entering flux from my neighbor will be distributed on my slip systems according to the !*compatibility considerEnteringFlux = .false. neighbor_v = 0.0_pReal ! needed for check of sign change in flux density below neighbor_rhoSgl = 0.0_pReal if (neighbor_n > 0_pInt) then if (phase_plasticity(material_phase(1,neighbor_ip,neighbor_el)) == PLASTICITY_NONLOCAL_ID & .and. any(compatibility(:,:,:,n,ip,el) > 0.0_pReal)) & considerEnteringFlux = .true. endif if (considerEnteringFlux) then forall (s = 1:ns, t = 1_pInt:4_pInt) 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 * mesh_ipVolume(neighbor_ip,neighbor_el) ** 0.667_pReal < prm%significantN & .or. neighbor_rhoSgl < prm%significantRho) & neighbor_rhoSgl = 0.0_pReal normal_neighbor2me_defConf = math_det33(Favg) * math_mul33x3(math_inv33(transpose(Favg)), & mesh_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 = math_mul33x3(transpose(neighbor_Fe), normal_neighbor2me_defConf) & / math_det33(neighbor_Fe) ! interface normal in the lattice configuration of my neighbor area = mesh_ipArea(neighbor_n,neighbor_ip,neighbor_el) * norm2(normal_neighbor2me) normal_neighbor2me = normal_neighbor2me / norm2(normal_neighbor2me) ! normalize the surface normal to unit length do s = 1_pInt,ns do t = 1_pInt,4_pInt c = (t + 1_pInt) / 2 topp = t + mod(t,2_pInt) - mod(t+1_pInt,2_pInt) if (neighbor_v(s,t) * math_mul3x3(m(1:3,s,t), normal_neighbor2me) > 0.0_pReal & ! flux from my neighbor to me == entering flux for me .and. v(s,t) * neighbor_v(s,t) >= 0.0_pReal ) then ! ... only if no sign change in flux density lineLength = neighbor_rhoSgl(s,t) * neighbor_v(s,t) & * math_mul3x3(m(1:3,s,t), normal_neighbor2me) * area ! positive line length that wants to enter through this interface where (compatibility(c,1_pInt:ns,s,n,ip,el) > 0.0_pReal) & ! positive compatibility... rhoDotFlux(1_pInt:ns,t) = rhoDotFlux(1_pInt:ns,t) & + lineLength / mesh_ipVolume(ip,el) & ! ... transferring to equally signed mobile dislocation type * compatibility(c,1_pInt:ns,s,n,ip,el) ** 2.0_pReal where (compatibility(c,1_pInt:ns,s,n,ip,el) < 0.0_pReal) & ! ..negative compatibility... rhoDotFlux(1_pInt:ns,topp) = rhoDotFlux(1_pInt:ns,topp) & + lineLength / mesh_ipVolume(ip,el) & ! ... transferring to opposite signed mobile dislocation type * compatibility(c,1_pInt:ns,s,n,ip,el) ** 2.0_pReal endif enddo enddo endif !* FLUX FROM ME TO MY NEIGHBOR !* This is not considered, if my opposite neighbor has a different constitutive law than nonlocal (still considered for nonlocal law with 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_pInt) then if (phase_plasticity(material_phase(1,opposite_ip,opposite_el)) /= PLASTICITY_NONLOCAL_ID) & considerLeavingFlux = .false. endif if (considerLeavingFlux) then !* timeSyncing mode: If the central ip has zero subfraction, always use "state0". This is needed in case of !* a synchronization step for the central ip, because then "state" contains the values at the end of the !* previously converged full time step. Also, if either me or my neighbor has zero subfraction, we have to !* use "state0" to make sure that fluxes on both sides of the (potential) timestep are equal. my_rhoSgl = rhoSgl my_v = v normal_me2neighbor_defConf = math_det33(Favg) & * math_mul33x3(math_inv33(transpose(Favg)), & mesh_ipAreaNormal(1:3,n,ip,el)) ! calculate the normal of the interface in (average) deformed configuration (pointing from me to my neighbor!!!) normal_me2neighbor = math_mul33x3(transpose(my_Fe), normal_me2neighbor_defConf) & / math_det33(my_Fe) ! interface normal in my lattice configuration area = mesh_ipArea(n,ip,el) * norm2(normal_me2neighbor) normal_me2neighbor = normal_me2neighbor / norm2(normal_me2neighbor) ! normalize the surface normal to unit length do s = 1_pInt,ns do t = 1_pInt,4_pInt c = (t + 1_pInt) / 2_pInt if (my_v(s,t) * math_mul3x3(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 transmissivity = sum(compatibility(c,1_pInt: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) & * math_mul3x3(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 / mesh_ipVolume(ip,el) ! subtract dislocation flux from current type rhoDotFlux(s,t+4_pInt) = rhoDotFlux(s,t+4_pInt) & + lineLength / mesh_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 enddo ! neighbor loop endif !**************************************************************************** !*** calculate dipole formation and annihilation !*** formation by glide do c = 1_pInt,2_pInt 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_pInt:2_pInt) & rhoDotAthermalAnnihilation(1:ns,c+8_pInt) = -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) & ! only fcc forall (s = 1:ns, prm%colinearSystem(s) > 0_pInt) & 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 ) & * prm%mu / ( 2.0_pReal * PI * (1.0_pReal-prm%nu) ) & * 2.0_pReal / ( dUpper(1:ns,1) + dLower(1:ns,1) ) forall (s = 1_pInt:ns, dUpper(s,1) > dLower(s,1)) & rhoDotThermalAnnihilation(s,9) = max(- 4.0_pReal * rhoDip(s,1) * vClimb(s) / (dUpper(s,1) - dLower(s,1)), & - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) & - rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have !**************************************************************************** !*** assign the rates of dislocation densities to my dotState !*** if evolution rates lead to negative densities, a cutback is enforced rhoDot = 0.0_pReal rhoDot = rhoDotFlux & + rhoDotMultiplication & + rhoDotSingle2DipoleGlide & + rhoDotAthermalAnnihilation & + rhoDotThermalAnnihilation results(instance)%rhoDotFlux(1:ns,1:8,o) = rhoDotFlux(1:ns,1:8) results(instance)%rhoDotMultiplication(1:ns,1:2,o) = rhoDotMultiplication(1:ns,[1,3]) results(instance)%rhoDotSingle2DipoleGlide(1:ns,1:2,o) = rhoDotSingle2DipoleGlide(1:ns,9:10) results(instance)%rhoDotAthermalAnnihilation(1:ns,1:2,o) = rhoDotAthermalAnnihilation(1:ns,9:10) results(instance)%rhoDotThermalAnnihilation(1:ns,1:2,o) = rhoDotThermalAnnihilation(1:ns,9:10) results(instance)%rhoDotEdgeJogs(1:ns,o) = 2.0_pReal * rhoDotThermalAnnihilation(1:ns,1) #ifdef DEBUG if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt & .and. ((debug_e == el .and. debug_i == ip)& .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then write(6,'(a,/,4(12x,12(e12.5,1x),/))') '<< CONST >> dislocation multiplication', & rhoDotMultiplication(1:ns,1:4) * timestep write(6,'(a,/,8(12x,12(e12.5,1x),/))') '<< CONST >> dislocation flux', & rhoDotFlux(1:ns,1:8) * timestep write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> dipole formation by glide', & rhoDotSingle2DipoleGlide * timestep write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> athermal dipole annihilation', & rhoDotAthermalAnnihilation * timestep write(6,'(a,/,2(12x,12(e12.5,1x),/))') '<< CONST >> thermally activated dipole annihilation', & rhoDotThermalAnnihilation(1:ns,9:10) * timestep write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> total density change', & rhoDot * timestep write(6,'(a,/,10(12x,12(f12.5,1x),/))') '<< CONST >> relative density change', & rhoDot(1:ns,1:8) * timestep / (abs(stt%rho(:,sgl))+1.0e-10), & rhoDot(1:ns,9:10) * timestep / (stt%rho(:,dip)+1.0e-10) write(6,*) endif #endif if ( any(stt%rho(:,mob) + rhoDot(1:ns,1:4) * timestep < -prm%aTolRho) & .or. any(stt%rho(:,dip) + rhoDot(1:ns,9:10) * timestep < -prm%aTolRho)) then #ifdef DEBUG if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt) then write(6,'(a,i5,a,i2)') '<< CONST >> evolution rate leads to negative density at el ',el,' ip ',ip write(6,'(a)') '<< CONST >> enforcing cutback !!!' endif #endif plasticState(p)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN) return else forall (s = 1:ns, t = 1_pInt:4_pInt) plasticState(p)%dotState(iRhoU(s,t,instance),o) = rhoDot(s,t) plasticState(p)%dotState(iRhoB(s,t,instance),o) = rhoDot(s,t+4_pInt) endforall forall (s = 1:ns, c = 1_pInt:2_pInt) & plasticState(p)%dotState(iRhoD(s,c,instance),o) = rhoDot(s,c+8_pInt) forall (s = 1:ns) & dot%accumulatedshear(s,o) = sum(gdot(s,1:4)) endif end associate end subroutine plastic_nonlocal_dotState !********************************************************************* !* COMPATIBILITY UPDATE * !* 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. * !********************************************************************* subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) use math, only: math_mul3x3, math_qRot use rotations, only: rotation use material, only: material_phase, & material_texture, & phase_localPlasticity, & phase_plasticityInstance use mesh, only: mesh_ipNeighborhood, & theMesh use lattice, only: lattice_qDisorientation implicit none !* input variables integer(pInt), intent(in) :: i, & ! ip index e ! element index type(rotation), dimension(1,theMesh%elem%nIPs,theMesh%nElems), intent(in) :: & orientation ! crystal orientation in quaternions !* local variables integer(pInt) Nneighbors, & ! number of neighbors n, & ! neighbor index neighbor_e, & ! element index of my neighbor neighbor_i, & ! integration point index of my neighbor ph, & neighbor_phase, & textureID, & neighbor_textureID, & instance, & ! instance of plasticity ns, & ! number of active slip systems s1, & ! slip system index (me) s2 ! slip system index (my neighbor) real(pReal), dimension(4) :: absoluteMisorientation ! absolute misorientation (without symmetry) between me and my neighbor real(pReal), dimension(2,totalNslip(phase_plasticityInstance(material_phase(1,i,e))),& totalNslip(phase_plasticityInstance(material_phase(1,i,e))),& theMesh%elem%nIPneighbors) :: & my_compatibility ! my_compatibility for current element and ip real(pReal) :: my_compatibilitySum, & thresholdValue, & nThresholdValues logical, dimension(totalNslip(phase_plasticityInstance(material_phase(1,i,e)))) :: & belowThreshold type(rotation) :: rot Nneighbors = theMesh%elem%nIPneighbors ph = material_phase(1,i,e) textureID = material_texture(1,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_pInt:ns) my_compatibility(1:2,s1,s1,1:Nneighbors) = 1.0_pReal !*** Loop thrugh neighbors and check whether there is any my_compatibility. neighbors: do n = 1_pInt,Nneighbors neighbor_e = mesh_ipNeighborhood(1,n,i,e) neighbor_i = mesh_ipNeighborhood(2,n,i,e) !* FREE SURFACE !* Set surface transmissivity to the value specified in the material.config if (neighbor_e <= 0_pInt .or. neighbor_i <= 0_pInt) then forall(s1 = 1_pInt:ns) my_compatibility(1:2,s1,s1,n) = sqrt(prm%surfaceTransmissivity) cycle endif !* PHASE BOUNDARY !* If we encounter a different nonlocal "cpfem" phase at the neighbor, !* we consider this to be a real "physical" phase boundary, so completely incompatible. !* If one of the two "CPFEM" phases has a local plasticity law, !* we do not consider this to be a phase boundary, so completely compatible. neighbor_phase = material_phase(1,neighbor_i,neighbor_e) if (neighbor_phase /= ph) then if (.not. phase_localPlasticity(neighbor_phase) .and. .not. phase_localPlasticity(ph))& forall(s1 = 1_pInt: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 neighbor_textureID = material_texture(1,neighbor_i,neighbor_e) if (neighbor_textureID /= textureID) then if (.not. phase_localPlasticity(neighbor_phase)) then forall(s1 = 1_pInt:ns) & my_compatibility(1:2,s1,s1,n) = sqrt(prm%grainboundaryTransmissivity) 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), !* only values above or equal to a certain threshold value are considered. This threshold value is chosen, such that !* the number of compatible slip systems is minimized with the sum of the original my_compatibility values exceeding one. !* Finally the smallest my_compatibility value is decreased until the sum is exactly equal to one. !* All values below the threshold are set to zero. else rot = orientation(1,i,e)%misorientation(orientation(1,neighbor_i,neighbor_e)) absoluteMisorientation = rot%asQuaternion() mySlipSystems: do s1 = 1_pInt,ns neighborSlipSystems: do s2 = 1_pInt,ns my_compatibility(1,s2,s1,n) = math_mul3x3(prm%slip_normal(1:3,s1), & math_qRot(absoluteMisorientation, prm%slip_normal(1:3,s2))) & * abs(math_mul3x3(prm%slip_direction(1:3,s1), & math_qRot(absoluteMisorientation, prm%slip_direction(1:3,s2)))) my_compatibility(2,s2,s1,n) = abs(math_mul3x3(prm%slip_normal(1:3,s1), & math_qRot(absoluteMisorientation, prm%slip_normal(1:3,s2)))) & * abs(math_mul3x3(prm%slip_direction(1:3,s1), & math_qRot(absoluteMisorientation, 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))) thresholdValue = maxval(my_compatibility(2,1:ns,s1,n), belowThreshold(1:ns)) ! screws always positive nThresholdValues = real(count(my_compatibility(2,1:ns,s1,n) >= thresholdValue),pReal) where (my_compatibility(2,1:ns,s1,n) >= thresholdValue) & belowThreshold(1:ns) = .false. if (my_compatibilitySum + thresholdValue * nThresholdValues > 1.0_pReal) & where (abs(my_compatibility(1:2,1:ns,s1,n)) >= thresholdValue) & ! MD: rather check below threshold? my_compatibility(1:2,1:ns,s1,n) = sign((1.0_pReal - my_compatibilitySum) & / nThresholdValues, my_compatibility(1:2,1:ns,s1,n)) my_compatibilitySum = my_compatibilitySum + nThresholdValues * thresholdValue enddo where (belowThreshold(1:ns)) my_compatibility(1,1:ns,s1,n) = 0.0_pReal 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 return array of constitutive results !-------------------------------------------------------------------------------------------------- function plastic_nonlocal_postResults(Mp,ip,el) result(postResults) use prec, only: & dNeq0 use math, only: & math_mul33x3, & math_mul33xx33, & pi use material, only: & material_phase, & phaseAt, phasememberAt, & plasticState, & phase_plasticityInstance implicit none real(pReal), dimension(3,3), intent(in) :: Mp !< MandelStress integer(pInt), intent(in) :: & ip, & !< integration point el !< element real(pReal), dimension(sum(plastic_nonlocal_sizePostResult(:,phase_plasticityInstance(material_phase(1_pInt,ip,el))))) :: & postResults integer(pInt) :: & ph, & instance, & !< current instance of this plasticity ns, & !< short notation for the total number of active slip systems c, & !< character of dislocation cs, & !< constitutive result index o, & !< index of current output of,& !< offset shortcut t, & !< type of dislocation s !< index of my current slip system real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1_pInt,ip,el))),8) :: & rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles) rhoDotSgl !< evolution rate of single dislocation densities (positive/negative screw and edge without dipoles) real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1_pInt,ip,el))),4) :: & gdot, & !< shear rates v !< velocities real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1_pInt,ip,el)))) :: & rhoForest, & !< forest dislocation density tau !< current resolved shear stress real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1_pInt,ip,el))),2) :: & rhoDip, & !< current dipole dislocation densities (screw and edge dipoles) rhoDotDip, & !< evolution rate of dipole dislocation densities (screw and edge dipoles) dLower, & !< minimum stable dipole distance for edges and screws dUpper !< current maximum stable dipole distance for edges and screws ph = phaseAt(1,ip,el) of = phasememberAt(1,ip,el) instance = phase_plasticityInstance(ph) ns = totalNslip(instance) cs = 0_pInt associate(prm => param(instance),dst => microstructure(instance),stt=>state(instance)) !* short hand notations for state variables forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) rhoSgl(s,t) = plasticState(ph)%State(iRhoU(s,t,instance),of) rhoSgl(s,t+4_pInt) = plasticState(ph)%State(iRhoB(s,t,instance),of) v(s,t) = plasticState(ph)%State(iV(s,t,instance),of) rhoDotSgl(s,t) = plasticState(ph)%dotState(iRhoU(s,t,instance),of) rhoDotSgl(s,t+4_pInt) = plasticState(ph)%dotState(iRhoB(s,t,instance),of) endforall forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) rhoDip(s,c) = plasticState(ph)%State(iRhoD(s,c,instance),of) rhoDotDip(s,c) = plasticState(ph)%dotState(iRhoD(s,c,instance),of) endforall rhoForest = plasticState(ph)%State(iRhoF(1:ns,instance),of) !* Calculate shear rate forall (t = 1_pInt:4_pInt) & gdot(1:ns,t) = rhoSgl(1:ns,t) * prm%burgers(1:ns) * v(1:ns,t) !* calculate limits for stable dipole height do s = 1_pInt,ns 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 dLower = prm%minDipoleHeight(1:ns,1:2) dUpper(1:ns,1) = prm%mu * prm%burgers(1:ns) & / (8.0_pReal * pi * (1.0_pReal - prm%nu) * abs(tau)) dUpper(1:ns,2) = prm%mu * prm%burgers(1:ns) & / (4.0_pReal * pi * abs(tau)) do c = 1, 2 where(dNeq0(sqrt(rhoSgl(1:ns,2*c-1)+rhoSgl(1:ns,2*c)+abs(rhoSgl(1:ns,2*c+3))& +abs(rhoSgl(1:ns,2*c+4))+rhoDip(1:ns,c)))) & dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) & + abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), & dUpper(1:ns,c)) enddo dUpper = max(dUpper,dLower) outputsLoop: do o = 1_pInt,size(param(instance)%outputID) select case(param(instance)%outputID(o)) case (rho_sgl_mob_edg_pos_ID) postResults(cs+1_pInt:cs+ns) = stt%rho_sgl_mob_edg_pos(:,of) cs = cs + ns case (rho_sgl_imm_edg_pos_ID) postResults(cs+1_pInt:cs+ns) = stt%rho_sgl_imm_edg_pos(:,of) cs = cs + ns case (rho_sgl_mob_edg_neg_ID) postResults(cs+1_pInt:cs+ns) = stt%rho_sgl_mob_edg_neg(:,of) cs = cs + ns case (rho_sgl_imm_edg_neg_ID) postResults(cs+1_pInt:cs+ns) = stt%rho_sgl_imm_edg_neg(:,of) cs = cs + ns case (rho_dip_edg_ID) postResults(cs+1_pInt:cs+ns) = stt%rho_dip_edg(:,of) cs = cs + ns case (rho_sgl_mob_scr_pos_ID) postResults(cs+1_pInt:cs+ns) = stt%rho_sgl_mob_scr_pos(:,of) cs = cs + ns case (rho_sgl_imm_scr_pos_ID) postResults(cs+1_pInt:cs+ns) = stt%rho_sgl_imm_scr_pos(:,of) cs = cs + ns case (rho_sgl_mob_scr_neg_ID) postResults(cs+1_pInt:cs+ns) = stt%rho_sgl_mob_scr_neg(:,of) cs = cs + ns case (rho_sgl_imm_scr_neg_ID) postResults(cs+1_pInt:cs+ns) = stt%rho_sgl_imm_scr_neg(:,of) cs = cs + ns case (rho_dip_scr_ID) postResults(cs+1_pInt:cs+ns) = stt%rho_dip_scr(:,of) cs = cs + ns case (rho_forest_ID) postResults(cs+1_pInt:cs+ns) = rhoForest cs = cs + ns case (shearrate_ID) postResults(cs+1_pInt:cs+ns) = sum(gdot,2) cs = cs + ns case (resolvedstress_ID) postResults(cs+1_pInt:cs+ns) = tau cs = cs + ns case (resolvedstress_back_ID) postResults(cs+1_pInt:cs+ns) = dst%tau_back(:,of) cs = cs + ns case (resolvedstress_external_ID) do s = 1_pInt,ns postResults(cs+s) = math_mul33xx33(Mp, prm%Schmid(1:3,1:3,s)) enddo cs = cs + ns case (resistance_ID) postResults(cs+1_pInt:cs+ns) = dst%tau_Threshold(:,of) cs = cs + ns case (rho_dot_sgl_ID) postResults(cs+1_pInt:cs+ns) = sum(rhoDotSgl(1:ns,1:4),2) & + sum(rhoDotSgl(1:ns,5:8)*sign(1.0_pReal,rhoSgl(1:ns,5:8)),2) cs = cs + ns case (rho_dot_sgl_mobile_ID) postResults(cs+1_pInt:cs+ns) = sum(rhoDotSgl(1:ns,1:4),2) cs = cs + ns case (rho_dot_dip_ID) postResults(cs+1_pInt:cs+ns) = sum(rhoDotDip,2) cs = cs + ns case (rho_dot_gen_ID) ! Obsolete postResults(cs+1_pInt:cs+ns) = results(instance)%rhoDotMultiplication(1:ns,1,of) & + results(instance)%rhoDotMultiplication(1:ns,2,of) cs = cs + ns case (rho_dot_gen_edge_ID) postResults(cs+1_pInt:cs+ns) = results(instance)%rhoDotMultiplication(1:ns,1,of) cs = cs + ns case (rho_dot_gen_screw_ID) postResults(cs+1_pInt:cs+ns) = results(instance)%rhoDotMultiplication(1:ns,2,of) cs = cs + ns case (rho_dot_sgl2dip_edge_ID) postResults(cs+1_pInt:cs+ns) = results(instance)%rhoDotSingle2DipoleGlide(1:ns,1,of) cs = cs + ns case (rho_dot_sgl2dip_screw_ID) postResults(cs+1_pInt:cs+ns) = results(instance)%rhoDotSingle2DipoleGlide(1:ns,2,of) cs = cs + ns case (rho_dot_ann_ath_ID) postResults(cs+1_pInt:cs+ns) = results(instance)%rhoDotAthermalAnnihilation(1:ns,1,of) & + results(instance)%rhoDotAthermalAnnihilation(1:ns,2,of) cs = cs + ns case (rho_dot_ann_the_edge_ID) postResults(cs+1_pInt:cs+ns) = results(instance)%rhoDotThermalAnnihilation(1:ns,1,of) cs = cs + ns case (rho_dot_ann_the_screw_ID) postResults(cs+1_pInt:cs+ns) = results(instance)%rhoDotThermalAnnihilation(1:ns,2,of) cs = cs + ns case (rho_dot_edgejogs_ID) postResults(cs+1_pInt:cs+ns) = results(instance)%rhoDotEdgeJogs(1:ns,of) cs = cs + ns case (rho_dot_flux_mobile_ID) postResults(cs+1_pInt:cs+ns) = sum(results(instance)%rhoDotFlux(1:ns,1:4,of),2) cs = cs + ns case (rho_dot_flux_edge_ID) postResults(cs+1_pInt:cs+ns) = sum(results(instance)%rhoDotFlux(1:ns,1:2,of),2) & + sum(results(instance)%rhoDotFlux(1:ns,5:6,of)*sign(1.0_pReal,rhoSgl(1:ns,5:6)),2) cs = cs + ns case (rho_dot_flux_screw_ID) postResults(cs+1_pInt:cs+ns) = sum(results(instance)%rhoDotFlux(1:ns,3:4,of),2) & + sum(results(instance)%rhoDotFlux(1:ns,7:8,of)*sign(1.0_pReal,rhoSgl(1:ns,7:8)),2) cs = cs + ns case (velocity_edge_pos_ID) postResults(cs+1_pInt:cs+ns) = v(1:ns,1) cs = cs + ns case (velocity_edge_neg_ID) postResults(cs+1_pInt:cs+ns) = v(1:ns,2) cs = cs + ns case (velocity_screw_pos_ID) postResults(cs+1_pInt:cs+ns) = v(1:ns,3) cs = cs + ns case (velocity_screw_neg_ID) postResults(cs+1_pInt:cs+ns) = v(1:ns,4) cs = cs + ns case (maximumdipoleheight_edge_ID) postResults(cs+1_pInt:cs+ns) = dUpper(1:ns,1) cs = cs + ns case (maximumdipoleheight_screw_ID) postResults(cs+1_pInt:cs+ns) = dUpper(1:ns,2) cs = cs + ns case(accumulatedshear_ID) postResults(cs+1_pInt:cs+ns) = stt%accumulatedshear(:,of) cs = cs + ns end select enddo outputsLoop end associate end function plastic_nonlocal_postResults function getRho(instance,of,ip,el) use mesh implicit none 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 mobile 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) * mesh_ipVolume(ip,el) ** 0.667_pReal < prm%significantN & .or. abs(getRho) < prm%significantRho) & getRho = 0.0_pReal end associate end function getRho !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- subroutine plastic_nonlocal_results(instance,group) #if defined(PETSc) || defined(DAMASKHDF5) use results implicit none integer, intent(in) :: instance character(len=*) :: group integer :: o associate(prm => param(instance), stt => state(instance)) outputsLoop: do o = 1_pInt,size(prm%outputID) select case(prm%outputID(o)) end select enddo outputsLoop end associate #else integer, intent(in) :: instance character(len=*) :: group #endif end subroutine plastic_nonlocal_results end module plastic_nonlocal