From fbc607e14d8f75f34b6178739ff7afc3974c4c2e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 18 Mar 2019 22:53:24 +0100 Subject: [PATCH 01/12] no pInt --- src/plastic_disloUCLA.f90 | 110 +++++++++++++++++++------------------- 1 file changed, 55 insertions(+), 55 deletions(-) diff --git a/src/plastic_disloUCLA.f90 b/src/plastic_disloUCLA.f90 index c67ff4ad0..817070b8d 100644 --- a/src/plastic_disloUCLA.f90 +++ b/src/plastic_disloUCLA.f90 @@ -12,9 +12,9 @@ module plastic_disloUCLA implicit none private - integer(pInt), dimension(:,:), allocatable, target, public :: & + integer, dimension(:,:), allocatable, target, public :: & plastic_disloUCLA_sizePostResult !< size of each post result output - character(len=64), dimension(:,:), allocatable, target, public :: & + character(len=64), dimension(:,:), allocatable, target, public :: & plastic_disloUCLA_output !< name of each post result output real(pReal), parameter, private :: & @@ -39,7 +39,7 @@ module plastic_disloUCLA mu, & D0, & !< prefactor for self-diffusion coefficient Qsd !< activation energy for dislocation climb - real(pReal), allocatable, dimension(:) :: & + real(pReal), dimension(:), allocatable :: & rho0, & !< initial edge dislocation density rhoDip0, & !< initial edge dipole density burgers, & !< absolute length of burgers vector [m] @@ -58,32 +58,32 @@ module plastic_disloUCLA kink_height, & !< height of the kink pair w, & !< width of the kink pair omega !< attempt frequency for kink pair nucleation - real(pReal), allocatable, dimension(:,:) :: & + real(pReal), dimension(:,:), allocatable :: & interaction_SlipSlip, & !< slip resistance from slip activity forestProjectionEdge - real(pReal), allocatable, dimension(:,:,:) :: & + real(pReal), dimension(:,:,:), allocatable :: & Schmid, & nonSchmid_pos, & nonSchmid_neg - integer(pInt) :: & + integer :: & totalNslip !< total number of active slip system - integer(pInt), allocatable, dimension(:) :: & + integer, dimension(:), allocatable, :: & Nslip !< number of active slip systems for each family - integer(kind(undefined_ID)), allocatable, dimension(:) :: & + integer(kind(undefined_ID)), dimension(:),allocatable :: & outputID !< ID of each post result output logical :: & dipoleFormation !< flag indicating consideration of dipole formation end type !< container type for internal constitutive parameters type, private :: tDisloUCLAState - real(pReal), pointer, dimension(:,:) :: & + real(pReal), dimension(:,:), pointer :: & rhoEdge, & rhoEdgeDip, & accshear end type tDisloUCLAState type, private :: tDisloUCLAdependentState - real(pReal), allocatable, dimension(:,:) :: & + real(pReal), dimension(:,:), allocatable :: & mfp, & dislocationSpacing, & threshold_stress @@ -139,14 +139,14 @@ subroutine plastic_disloUCLA_init() use lattice implicit none - integer(pInt) :: & + integer :: & Ninstance, & p, i, & NipcMyPhase, & sizeState, sizeDotState, & startIndex, endIndex - integer(pInt), dimension(0), parameter :: emptyIntArray = [integer(pInt)::] + integer, dimension(0), parameter :: emptyIntArray = [integer::] real(pReal), dimension(0), parameter :: emptyRealArray = [real(pReal)::] character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] @@ -164,10 +164,10 @@ subroutine plastic_disloUCLA_init() write(6,'(a)') ' https://dx.doi.org/10.1016/j.ijplas.2015.09.002' Ninstance = int(count(phase_plasticity == PLASTICITY_DISLOUCLA_ID),pInt) - if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & + if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - allocate(plastic_disloUCLA_sizePostResult(maxval(phase_Noutput),Ninstance),source=0_pInt) + allocate(plastic_disloUCLA_sizePostResult(maxval(phase_Noutput),Ninstance),source=0) allocate(plastic_disloUCLA_output(maxval(phase_Noutput),Ninstance)) plastic_disloUCLA_output = '' @@ -177,7 +177,7 @@ subroutine plastic_disloUCLA_init() allocate(dependentState(Ninstance)) - do p = 1_pInt, size(phase_plasticity) + do p = 1, size(phase_plasticity) if (phase_plasticity(p) /= PLASTICITY_DISLOUCLA_ID) cycle associate(prm => param(phase_plasticityInstance(p)), & dot => dotState(phase_plasticityInstance(p)), & @@ -198,15 +198,15 @@ subroutine plastic_disloUCLA_init() ! slip related parameters prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) prm%totalNslip = sum(prm%Nslip) - slipActive: if (prm%totalNslip > 0_pInt) then + slipActive: if (prm%totalNslip > 0) 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) + prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1) + prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1) else prm%nonSchmid_pos = prm%Schmid prm%nonSchmid_neg = prm%Schmid @@ -227,9 +227,9 @@ subroutine plastic_disloUCLA_init() prm%clambda = config%getFloats('clambdaslip', requiredSize=size(prm%Nslip)) prm%tau_Peierls = config%getFloats('tau_peierls', requiredSize=size(prm%Nslip)) ! ToDo: Deprecated prm%p = config%getFloats('p_slip', requiredSize=size(prm%Nslip), & - defaultVal=[(1.0_pReal,i=1_pInt,size(prm%Nslip))]) + defaultVal=[(1.0_pReal,i=1,size(prm%Nslip))]) prm%q = config%getFloats('q_slip', requiredSize=size(prm%Nslip), & - defaultVal=[(1.0_pReal,i=1_pInt,size(prm%Nslip))]) + defaultVal=[(1.0_pReal,i=1,size(prm%Nslip))]) prm%kink_height = config%getFloats('kink_height', requiredSize=size(prm%Nslip)) prm%w = config%getFloats('kink_width', requiredSize=size(prm%Nslip)) prm%omega = config%getFloats('omega', requiredSize=size(prm%Nslip)) @@ -282,28 +282,28 @@ subroutine plastic_disloUCLA_init() !-------------------------------------------------------------------------------------------------- ! exit if any parameter is out of range if (extmsg /= '') & - call IO_error(211_pInt,ext_msg=trim(extmsg)//'('//PLASTICITY_DISLOUCLA_label//')') + call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_DISLOUCLA_label//')') !-------------------------------------------------------------------------------------------------- ! output pararameters outputs = config%getStrings('(output)',defaultVal=emptyStringArray) allocate(prm%outputID(0)) - do i=1_pInt, size(outputs) + do i=1, size(outputs) outputID = undefined_ID select case(trim(outputs(i))) case ('edge_density') - outputID = merge(rho_ID,undefined_ID,prm%totalNslip>0_pInt) + outputID = merge(rho_ID,undefined_ID,prm%totalNslip>0) case ('dipole_density') - outputID = merge(rhoDip_ID,undefined_ID,prm%totalNslip>0_pInt) + outputID = merge(rhoDip_ID,undefined_ID,prm%totalNslip>0) case ('shear_rate','shearrate','shear_rate_slip','shearrate_slip') - outputID = merge(shearrate_ID,undefined_ID,prm%totalNslip>0_pInt) + outputID = merge(shearrate_ID,undefined_ID,prm%totalNslip>0) case ('accumulated_shear','accumulatedshear','accumulated_shear_slip') - outputID = merge(accumulatedshear_ID,undefined_ID,prm%totalNslip>0_pInt) + outputID = merge(accumulatedshear_ID,undefined_ID,prm%totalNslip>0) case ('mfp','mfp_slip') - outputID = merge(mfp_ID,undefined_ID,prm%totalNslip>0_pInt) + outputID = merge(mfp_ID,undefined_ID,prm%totalNslip>0) case ('threshold_stress','threshold_stress_slip') - outputID = merge(thresholdstress_ID,undefined_ID,prm%totalNslip>0_pInt) + outputID = merge(thresholdstress_ID,undefined_ID,prm%totalNslip>0) end select @@ -318,30 +318,30 @@ subroutine plastic_disloUCLA_init() !-------------------------------------------------------------------------------------------------- ! allocate state arrays NipcMyPhase = count(material_phase == p) - sizeDotState = int(size(['rhoEdge ','rhoEdgeDip ','accshearslip']),pInt) * prm%totalNslip + sizeDotState = size(['rhoEdge ','rhoEdgeDip ','accshearslip']) * prm%totalNslip sizeState = sizeDotState - call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0_pInt, & - prm%totalNslip,0_pInt,0_pInt) + call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0, & + prm%totalNslip,0,0) plasticState(p)%sizePostResults = sum(plastic_disloUCLA_sizePostResult(:,phase_plasticityInstance(p))) !-------------------------------------------------------------------------------------------------- ! locally defined state aliases and initialization of state0 and aTolState - startIndex = 1_pInt + startIndex = 1 endIndex = prm%totalNslip stt%rhoEdge=>plasticState(p)%state(startIndex:endIndex,:) stt%rhoEdge= spread(prm%rho0,2,NipcMyPhase) dot%rhoEdge=>plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolRho - startIndex = endIndex + 1_pInt + startIndex = endIndex + 1 endIndex = endIndex + prm%totalNslip stt%rhoEdgeDip=>plasticState(p)%state(startIndex:endIndex,:) stt%rhoEdgeDip= spread(prm%rhoDip0,2,NipcMyPhase) dot%rhoEdgeDip=>plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolRho - startIndex = endIndex + 1_pInt + startIndex = endIndex + 1 endIndex = endIndex + prm%totalNslip stt%accshear=>plasticState(p)%state(startIndex:endIndex,:) dot%accshear=>plasticState(p)%dotState(startIndex:endIndex,:) @@ -378,11 +378,11 @@ pure subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp,Mp,Temperature,inst Mp !< Mandel stress real(pReal), intent(in) :: & temperature !< temperature - integer(pInt), intent(in) :: & + integer, intent(in) :: & instance, & of - integer(pInt) :: & + integer :: & i,k,l,m,n real(pReal), dimension(param(instance)%totalNslip) :: & gdot_pos,gdot_neg, & @@ -394,9 +394,9 @@ pure subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp,Mp,Temperature,inst associate(prm => param(instance)) call kinetics(Mp,Temperature,instance,of,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) - do i = 1_pInt, prm%totalNslip + do i = 1, prm%totalNslip Lp = Lp + (gdot_pos(i)+gdot_neg(i))*prm%Schmid(1:3,1:3,i) - forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & + forall (k=1:3,l=1:3,m=1:3,n=1:3) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + dgdot_dtau_pos(i) * prm%Schmid(k,l,i) * prm%nonSchmid_pos(m,n,i) & + dgdot_dtau_neg(i) * prm%Schmid(k,l,i) * prm%nonSchmid_neg(m,n,i) @@ -423,7 +423,7 @@ subroutine plastic_disloUCLA_dotState(Mp,Temperature,instance,of) Mp !< Mandel stress real(pReal), intent(in) :: & temperature !< temperature - integer(pInt), intent(in) :: & + integer, intent(in) :: & instance, & of @@ -478,16 +478,16 @@ end subroutine plastic_disloUCLA_dotState subroutine plastic_disloUCLA_dependentState(instance,of) implicit none - integer(pInt), intent(in) :: & + integer, intent(in) :: & instance, & of - integer(pInt) :: & + integer :: & i associate(prm => param(instance), stt => state(instance),dst => dependentState(instance)) - forall (i = 1_pInt:prm%totalNslip) + forall (i = 1:prm%totalNslip) dst%dislocationSpacing(i,of) = sqrt(dot_product(stt%rhoEdge(:,of)+stt%rhoEdgeDip(:,of), & prm%forestProjectionEdge(:,i))) dst%threshold_stress(i,of) = prm%mu*prm%burgers(i) & @@ -518,38 +518,38 @@ function plastic_disloUCLA_postResults(Mp,Temperature,instance,of) result(postRe Mp !< Mandel stress real(pReal), intent(in) :: & temperature !< temperature - integer(pInt), intent(in) :: & + integer, intent(in) :: & instance, & of real(pReal), dimension(sum(plastic_disloUCLA_sizePostResult(:,instance))) :: & postResults - integer(pInt) :: & + integer :: & o,c,i real(pReal), dimension(param(instance)%totalNslip) :: & gdot_pos,gdot_neg - c = 0_pInt + c = 0 associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) - outputsLoop: do o = 1_pInt,size(prm%outputID) + outputsLoop: do o = 1,size(prm%outputID) select case(prm%outputID(o)) case (rho_ID) - postResults(c+1_pInt:c+prm%totalNslip) = stt%rhoEdge(1_pInt:prm%totalNslip,of) + postResults(c+1:c+prm%totalNslip) = stt%rhoEdge(1:prm%totalNslip,of) case (rhoDip_ID) - postResults(c+1_pInt:c+prm%totalNslip) = stt%rhoEdgeDip(1_pInt:prm%totalNslip,of) + postResults(c+1:c+prm%totalNslip) = stt%rhoEdgeDip(1:prm%totalNslip,of) case (shearrate_ID) call kinetics(Mp,Temperature,instance,of,gdot_pos,gdot_neg) postResults(c+1:c+prm%totalNslip) = gdot_pos + gdot_neg case (accumulatedshear_ID) - postResults(c+1_pInt:c+prm%totalNslip) = stt%accshear(1_pInt:prm%totalNslip, of) + postResults(c+1:c+prm%totalNslip) = stt%accshear(1:prm%totalNslip, of) case (mfp_ID) - postResults(c+1_pInt:c+prm%totalNslip) = dst%mfp(1_pInt:prm%totalNslip, of) + postResults(c+1:c+prm%totalNslip) = dst%mfp(1:prm%totalNslip, of) case (thresholdstress_ID) - postResults(c+1_pInt:c+prm%totalNslip) = dst%threshold_stress(1_pInt:prm%totalNslip,of) + postResults(c+1:c+prm%totalNslip) = dst%threshold_stress(1:prm%totalNslip,of) end select @@ -575,7 +575,7 @@ subroutine plastic_disloUCLA_results(instance,group) integer :: o associate(prm => param(instance), stt => state(instance)) - outputsLoop: do o = 1_pInt,size(prm%outputID) + outputsLoop: do o = 1,size(prm%outputID) select case(prm%outputID(o)) end select enddo outputsLoop @@ -609,7 +609,7 @@ pure subroutine kinetics(Mp,Temperature,instance,of, & Mp !< Mandel stress real(pReal), intent(in) :: & temperature !< temperature - integer(pInt), intent(in) :: & + integer, intent(in) :: & instance, & of @@ -627,11 +627,11 @@ pure subroutine kinetics(Mp,Temperature,instance,of, & dvel, vel, & tau_pos,tau_neg, & needsGoodName ! ToDo: @Karo: any idea? - integer(pInt) :: j + integer :: j associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) - do j = 1_pInt, prm%totalNslip + do j = 1, prm%totalNslip tau_pos(j) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,j)) tau_neg(j) = math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,j)) enddo From f9e47a94aa1993dd4c35c71d50410f1e4e54d874 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 18 Mar 2019 22:57:23 +0100 Subject: [PATCH 02/12] no pInt --- src/plastic_disloUCLA.f90 | 39 +++++++++++++++++++-------------------- 1 file changed, 19 insertions(+), 20 deletions(-) diff --git a/src/plastic_disloUCLA.f90 b/src/plastic_disloUCLA.f90 index 817070b8d..99a1eb8a8 100644 --- a/src/plastic_disloUCLA.f90 +++ b/src/plastic_disloUCLA.f90 @@ -7,14 +7,13 @@ !-------------------------------------------------------------------------------------------------- module plastic_disloUCLA use prec, only: & - pReal, & - pInt + pReal implicit none private integer, dimension(:,:), allocatable, target, public :: & plastic_disloUCLA_sizePostResult !< size of each post result output - character(len=64), dimension(:,:), allocatable, target, public :: & + character(len=64), dimension(:,:), allocatable, target, public :: & plastic_disloUCLA_output !< name of each post result output real(pReal), parameter, private :: & @@ -61,13 +60,13 @@ module plastic_disloUCLA real(pReal), dimension(:,:), allocatable :: & interaction_SlipSlip, & !< slip resistance from slip activity forestProjectionEdge - real(pReal), dimension(:,:,:), allocatable :: & + real(pReal), dimension(:,:,:), allocatable :: & Schmid, & nonSchmid_pos, & nonSchmid_neg integer :: & totalNslip !< total number of active slip system - integer, dimension(:), allocatable, :: & + integer, dimension(:), allocatable, :: & Nslip !< number of active slip systems for each family integer(kind(undefined_ID)), dimension(:),allocatable :: & outputID !< ID of each post result output @@ -79,7 +78,7 @@ module plastic_disloUCLA real(pReal), dimension(:,:), pointer :: & rhoEdge, & rhoEdgeDip, & - accshear + gamma_sl end type tDisloUCLAState type, private :: tDisloUCLAdependentState @@ -146,9 +145,9 @@ subroutine plastic_disloUCLA_init() sizeState, sizeDotState, & startIndex, endIndex - integer, dimension(0), parameter :: emptyIntArray = [integer::] - real(pReal), dimension(0), parameter :: emptyRealArray = [real(pReal)::] - character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] + integer, dimension(0), parameter :: emptyIntArray = [integer::] + real(pReal), dimension(0), parameter :: emptyRealArray = [real(pReal)::] + character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] integer(kind(undefined_ID)) :: & outputID @@ -163,7 +162,7 @@ subroutine plastic_disloUCLA_init() write(6,'(/,a)') ' Cereceda et al., International Journal of Plasticity 78:242–256, 2016' write(6,'(a)') ' https://dx.doi.org/10.1016/j.ijplas.2015.09.002' - Ninstance = int(count(phase_plasticity == PLASTICITY_DISLOUCLA_ID),pInt) + Ninstance = count(phase_plasticity == PLASTICITY_DISLOUCLA_ID) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance @@ -343,8 +342,8 @@ subroutine plastic_disloUCLA_init() startIndex = endIndex + 1 endIndex = endIndex + prm%totalNslip - stt%accshear=>plasticState(p)%state(startIndex:endIndex,:) - dot%accshear=>plasticState(p)%dotState(startIndex:endIndex,:) + stt%gamma_sl=>plasticState(p)%state(startIndex:endIndex,:) + dot%gamma_sl=>plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = 1.0e6_pReal !ToDo: better make optional parameter ! global alias plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) @@ -376,9 +375,9 @@ pure subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp,Mp,Temperature,inst real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress - real(pReal), intent(in) :: & + real(pReal), intent(in) :: & temperature !< temperature - integer, intent(in) :: & + integer, intent(in) :: & instance, & of @@ -442,7 +441,7 @@ subroutine plastic_disloUCLA_dotState(Mp,Temperature,instance,of) gdot_pos,gdot_neg, & tau_pos1 = tau_pos,tau_neg1 = tau_neg) - dot%accshear(:,of) = (gdot_pos+gdot_neg) ! ToDo: needs to be abs + dot%gamma_sl(:,of) = (gdot_pos+gdot_neg) ! ToDo: needs to be abs VacancyDiffusion = prm%D0*exp(-prm%Qsd/(kB*Temperature)) where(dEq0(tau_pos)) ! ToDo: use avg of pos and neg @@ -452,7 +451,7 @@ subroutine plastic_disloUCLA_dotState(Mp,Temperature,instance,of) EdgeDipDistance = math_clip((3.0_pReal*prm%mu*prm%burgers)/(16.0_pReal*PI*abs(tau_pos)), & prm%minDipDistance, & ! lower limit dst%mfp(:,of)) ! upper limit - DotRhoDipFormation = merge(((2.0_pReal*EdgeDipDistance)/prm%burgers)* stt%rhoEdge(:,of)*abs(dot%accshear(:,of)), & ! ToDo: ignore region of spontaneous annihilation + DotRhoDipFormation = merge(((2.0_pReal*EdgeDipDistance)/prm%burgers)* stt%rhoEdge(:,of)*abs(dot%gamma_sl(:,of)), & ! ToDo: ignore region of spontaneous annihilation 0.0_pReal, & prm%dipoleformation) ClimbVelocity = (3.0_pReal*prm%mu*VacancyDiffusion*prm%atomicVolume/(2.0_pReal*pi*kB*Temperature)) & @@ -460,11 +459,11 @@ subroutine plastic_disloUCLA_dotState(Mp,Temperature,instance,of) DotRhoEdgeDipClimb = (4.0_pReal*ClimbVelocity*stt%rhoEdgeDip(:,of))/(EdgeDipDistance-prm%minDipDistance) ! ToDo: Discuss with Franz: Stress dependency? end where - dot%rhoEdge(:,of) = abs(dot%accshear(:,of))/(prm%burgers*dst%mfp(:,of)) & ! multiplication + dot%rhoEdge(:,of) = abs(dot%gamma_sl(:,of))/(prm%burgers*dst%mfp(:,of)) & ! multiplication - DotRhoDipFormation & - - (2.0_pReal*prm%minDipDistance)/prm%burgers*stt%rhoEdge(:,of)*abs(dot%accshear(:,of)) !* Spontaneous annihilation of 2 single edge dislocations + - (2.0_pReal*prm%minDipDistance)/prm%burgers*stt%rhoEdge(:,of)*abs(dot%gamma_sl(:,of)) !* Spontaneous annihilation of 2 single edge dislocations dot%rhoEdgeDip(:,of) = DotRhoDipFormation & - - (2.0_pReal*prm%minDipDistance)/prm%burgers*stt%rhoEdgeDip(:,of)*abs(dot%accshear(:,of)) & !* Spontaneous annihilation of a single edge dislocation with a dipole constituent + - (2.0_pReal*prm%minDipDistance)/prm%burgers*stt%rhoEdgeDip(:,of)*abs(dot%gamma_sl(:,of)) & !* Spontaneous annihilation of a single edge dislocation with a dipole constituent - DotRhoEdgeDipClimb end associate @@ -545,7 +544,7 @@ function plastic_disloUCLA_postResults(Mp,Temperature,instance,of) result(postRe call kinetics(Mp,Temperature,instance,of,gdot_pos,gdot_neg) postResults(c+1:c+prm%totalNslip) = gdot_pos + gdot_neg case (accumulatedshear_ID) - postResults(c+1:c+prm%totalNslip) = stt%accshear(1:prm%totalNslip, of) + postResults(c+1:c+prm%totalNslip) = stt%gamma_sl(1:prm%totalNslip, of) case (mfp_ID) postResults(c+1:c+prm%totalNslip) = dst%mfp(1:prm%totalNslip, of) case (thresholdstress_ID) From fc997554baf4bda9cd83dbe5b44c21cdbefb8e52 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 18 Mar 2019 23:11:28 +0100 Subject: [PATCH 03/12] Simplified tangent calculation, thx to Karo --- PRIVATE | 2 +- src/plastic_disloUCLA.f90 | 143 ++++++++++++++++---------------------- 2 files changed, 60 insertions(+), 85 deletions(-) diff --git a/PRIVATE b/PRIVATE index c79dc5f1b..5272b99b7 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit c79dc5f1be75f90b0638c230d56c962bfd3b2474 +Subproject commit 5272b99b734600fbd1538b32fbcbbbb5e3dbcd30 diff --git a/src/plastic_disloUCLA.f90 b/src/plastic_disloUCLA.f90 index 99a1eb8a8..b099518ab 100644 --- a/src/plastic_disloUCLA.f90 +++ b/src/plastic_disloUCLA.f90 @@ -39,7 +39,7 @@ module plastic_disloUCLA D0, & !< prefactor for self-diffusion coefficient Qsd !< activation energy for dislocation climb real(pReal), dimension(:), allocatable :: & - rho0, & !< initial edge dislocation density + rho_mob_0, & !< initial edge dislocation density rhoDip0, & !< initial edge dipole density burgers, & !< absolute length of burgers vector [m] nonSchmidCoeff, & @@ -66,7 +66,7 @@ module plastic_disloUCLA nonSchmid_neg integer :: & totalNslip !< total number of active slip system - integer, dimension(:), allocatable, :: & + integer, dimension(:), allocatable :: & Nslip !< number of active slip systems for each family integer(kind(undefined_ID)), dimension(:),allocatable :: & outputID !< ID of each post result output @@ -76,8 +76,8 @@ module plastic_disloUCLA type, private :: tDisloUCLAState real(pReal), dimension(:,:), pointer :: & - rhoEdge, & - rhoEdgeDip, & + rho_mob, & + rho_dip, & gamma_sl end type tDisloUCLAState @@ -217,7 +217,7 @@ subroutine plastic_disloUCLA_init() prm%forestProjectionEdge = lattice_forestProjection(prm%Nslip,config%getString('lattice_structure'),& config%getFloat('c/a',defaultVal=0.0_pReal)) - prm%rho0 = config%getFloats('rhoedge0', requiredSize=size(prm%Nslip)) + prm%rho_mob_0 = config%getFloats('rhoedge0', requiredSize=size(prm%Nslip)) prm%rhoDip0 = config%getFloats('rhoedgedip0', requiredSize=size(prm%Nslip)) prm%v0 = config%getFloats('v0', requiredSize=size(prm%Nslip)) prm%burgers = config%getFloats('slipburgers', requiredSize=size(prm%Nslip)) @@ -243,7 +243,7 @@ subroutine plastic_disloUCLA_init() prm%dipoleformation = config%getFloat('dipoleformationfactor') > 0.0_pReal !should be on by default, ToDo: change to /key/-type key ! expand: family => system - prm%rho0 = math_expand(prm%rho0, prm%Nslip) + prm%rho_mob_0 = math_expand(prm%rho_mob_0, prm%Nslip) prm%rhoDip0 = math_expand(prm%rhoDip0, prm%Nslip) prm%q = math_expand(prm%q, prm%Nslip) prm%p = math_expand(prm%p, prm%Nslip) @@ -264,7 +264,7 @@ subroutine plastic_disloUCLA_init() ! sanity checks if ( prm%D0 <= 0.0_pReal) extmsg = trim(extmsg)//' d0' if ( prm%Qsd <= 0.0_pReal) extmsg = trim(extmsg)//' qsd' - if (any(prm%rho0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoedge0' + if (any(prm%rho_mob_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoedge0' if (any(prm%rhoDip0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoedgedip0' if (any(prm%v0 < 0.0_pReal)) extmsg = trim(extmsg)//' v0' if (any(prm%burgers <= 0.0_pReal)) extmsg = trim(extmsg)//' slipburgers' @@ -274,7 +274,7 @@ subroutine plastic_disloUCLA_init() if (any(prm%atomicVolume <= 0.0_pReal)) extmsg = trim(extmsg)//' catomicvolume or slipburgers' else slipActive - allocate(prm%rho0(0)) + allocate(prm%rho_mob_0(0)) allocate(prm%rhoDip0(0)) endif slipActive @@ -328,16 +328,16 @@ subroutine plastic_disloUCLA_init() ! locally defined state aliases and initialization of state0 and aTolState startIndex = 1 endIndex = prm%totalNslip - stt%rhoEdge=>plasticState(p)%state(startIndex:endIndex,:) - stt%rhoEdge= spread(prm%rho0,2,NipcMyPhase) - dot%rhoEdge=>plasticState(p)%dotState(startIndex:endIndex,:) + stt%rho_mob=>plasticState(p)%state(startIndex:endIndex,:) + stt%rho_mob= spread(prm%rho_mob_0,2,NipcMyPhase) + dot%rho_mob=>plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolRho startIndex = endIndex + 1 endIndex = endIndex + prm%totalNslip - stt%rhoEdgeDip=>plasticState(p)%state(startIndex:endIndex,:) - stt%rhoEdgeDip= spread(prm%rhoDip0,2,NipcMyPhase) - dot%rhoEdgeDip=>plasticState(p)%dotState(startIndex:endIndex,:) + stt%rho_dip=>plasticState(p)%state(startIndex:endIndex,:) + stt%rho_dip= spread(prm%rhoDip0,2,NipcMyPhase) + dot%rho_dip=>plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolRho startIndex = endIndex + 1 @@ -451,19 +451,19 @@ subroutine plastic_disloUCLA_dotState(Mp,Temperature,instance,of) EdgeDipDistance = math_clip((3.0_pReal*prm%mu*prm%burgers)/(16.0_pReal*PI*abs(tau_pos)), & prm%minDipDistance, & ! lower limit dst%mfp(:,of)) ! upper limit - DotRhoDipFormation = merge(((2.0_pReal*EdgeDipDistance)/prm%burgers)* stt%rhoEdge(:,of)*abs(dot%gamma_sl(:,of)), & ! ToDo: ignore region of spontaneous annihilation + DotRhoDipFormation = merge(((2.0_pReal*EdgeDipDistance)/prm%burgers)* stt%rho_mob(:,of)*abs(dot%gamma_sl(:,of)), & ! ToDo: ignore region of spontaneous annihilation 0.0_pReal, & prm%dipoleformation) ClimbVelocity = (3.0_pReal*prm%mu*VacancyDiffusion*prm%atomicVolume/(2.0_pReal*pi*kB*Temperature)) & * (1.0_pReal/(EdgeDipDistance+prm%minDipDistance)) - DotRhoEdgeDipClimb = (4.0_pReal*ClimbVelocity*stt%rhoEdgeDip(:,of))/(EdgeDipDistance-prm%minDipDistance) ! ToDo: Discuss with Franz: Stress dependency? + DotRhoEdgeDipClimb = (4.0_pReal*ClimbVelocity*stt%rho_dip(:,of))/(EdgeDipDistance-prm%minDipDistance) ! ToDo: Discuss with Franz: Stress dependency? end where - dot%rhoEdge(:,of) = abs(dot%gamma_sl(:,of))/(prm%burgers*dst%mfp(:,of)) & ! multiplication + dot%rho_mob(:,of) = abs(dot%gamma_sl(:,of))/(prm%burgers*dst%mfp(:,of)) & ! multiplication - DotRhoDipFormation & - - (2.0_pReal*prm%minDipDistance)/prm%burgers*stt%rhoEdge(:,of)*abs(dot%gamma_sl(:,of)) !* Spontaneous annihilation of 2 single edge dislocations - dot%rhoEdgeDip(:,of) = DotRhoDipFormation & - - (2.0_pReal*prm%minDipDistance)/prm%burgers*stt%rhoEdgeDip(:,of)*abs(dot%gamma_sl(:,of)) & !* Spontaneous annihilation of a single edge dislocation with a dipole constituent + - (2.0_pReal*prm%minDipDistance)/prm%burgers*stt%rho_mob(:,of)*abs(dot%gamma_sl(:,of)) !* Spontaneous annihilation of 2 single edge dislocations + dot%rho_dip(:,of) = DotRhoDipFormation & + - (2.0_pReal*prm%minDipDistance)/prm%burgers*stt%rho_dip(:,of)*abs(dot%gamma_sl(:,of)) & !* Spontaneous annihilation of a single edge dislocation with a dipole constituent - DotRhoEdgeDipClimb end associate @@ -487,10 +487,10 @@ subroutine plastic_disloUCLA_dependentState(instance,of) associate(prm => param(instance), stt => state(instance),dst => dependentState(instance)) forall (i = 1:prm%totalNslip) - dst%dislocationSpacing(i,of) = sqrt(dot_product(stt%rhoEdge(:,of)+stt%rhoEdgeDip(:,of), & + dst%dislocationSpacing(i,of) = sqrt(dot_product(stt%rho_mob(:,of)+stt%rho_dip(:,of), & prm%forestProjectionEdge(:,i))) dst%threshold_stress(i,of) = prm%mu*prm%burgers(i) & - * sqrt(dot_product(stt%rhoEdge(:,of)+stt%rhoEdgeDip(:,of), & + * sqrt(dot_product(stt%rho_mob(:,of)+stt%rho_dip(:,of), & prm%interaction_SlipSlip(:,i))) end forall @@ -537,9 +537,9 @@ function plastic_disloUCLA_postResults(Mp,Temperature,instance,of) result(postRe select case(prm%outputID(o)) case (rho_ID) - postResults(c+1:c+prm%totalNslip) = stt%rhoEdge(1:prm%totalNslip,of) + postResults(c+1:c+prm%totalNslip) = stt%rho_mob(1:prm%totalNslip,of) case (rhoDip_ID) - postResults(c+1:c+prm%totalNslip) = stt%rhoEdgeDip(1:prm%totalNslip,of) + postResults(c+1:c+prm%totalNslip) = stt%rho_dip(1:prm%totalNslip,of) case (shearrate_ID) call kinetics(Mp,Temperature,instance,of,gdot_pos,gdot_neg) postResults(c+1:c+prm%totalNslip) = gdot_pos + gdot_neg @@ -595,7 +595,7 @@ end subroutine plastic_disloUCLA_results ! have the optional arguments at the end !-------------------------------------------------------------------------------------------------- pure subroutine kinetics(Mp,Temperature,instance,of, & - gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg,tau_pos1,tau_neg1) + gamma_pos,gamma_neg,dgamma_dtau_pos,dgamma_dtau_neg,tau_pos1,tau_neg1) use prec, only: & tol_math_check, & dEq, dNeq0 @@ -613,11 +613,11 @@ pure subroutine kinetics(Mp,Temperature,instance,of, & of real(pReal), intent(out), dimension(param(instance)%totalNslip) :: & - gdot_pos, & - gdot_neg + gamma_pos, & + gamma_neg real(pReal), intent(out), optional, dimension(param(instance)%totalNslip) :: & - dgdot_dtau_pos, & - dgdot_dtau_neg, & + dgamma_dtau_pos, & + dgamma_dtau_neg, & tau_pos1, & tau_neg1 real(pReal), dimension(param(instance)%totalNslip) :: & @@ -625,6 +625,7 @@ pure subroutine kinetics(Mp,Temperature,instance,of, & StressRatio_p,StressRatio_pminus1, & dvel, vel, & tau_pos,tau_neg, & + t_n, t_k, dtk,dtn, & needsGoodName ! ToDo: @Karo: any idea? integer :: j @@ -640,7 +641,7 @@ pure subroutine kinetics(Mp,Temperature,instance,of, & if (present(tau_neg1)) tau_neg1 = tau_neg associate(BoltzmannRatio => prm%H0kp/(kB*Temperature), & - DotGamma0 => stt%rhoEdge(:,of)*prm%burgers*prm%v0, & + DotGamma0 => stt%rho_mob(:,of)*prm%burgers*prm%v0, & effectiveLength => dst%mfp(:,of) - prm%w) significantPositiveTau: where(abs(tau_pos)-dst%threshold_stress(:,of) > tol_math_check) @@ -648,41 +649,28 @@ pure subroutine kinetics(Mp,Temperature,instance,of, & StressRatio_p = StressRatio** prm%p StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) + + t_n = prm%burgers/(needsGoodName*prm%omega*effectiveLength) + t_k = effectiveLength * prm%B /(2.0_pReal*prm%burgers*tau_pos) ! our definition of tk is different with the one in dislotwin - vel = 2.0_pReal*prm%burgers * prm%kink_height * prm%omega & - * effectiveLength * tau_pos * needsGoodName & - / ( 2.0_pReal*(prm%burgers**2.0_pReal)*tau_pos & - + prm%omega * prm%B * effectiveLength**2.0_pReal* needsGoodName & - ) + vel = prm%kink_height/(t_n + t_k) - gdot_pos = DotGamma0 * sign(vel,tau_pos) * 0.5_pReal + gamma_pos = DotGamma0 * sign(vel,tau_pos) * 0.5_pReal else where significantPositiveTau - gdot_pos = 0.0_pReal + gamma_pos = 0.0_pReal end where significantPositiveTau - if (present(dgdot_dtau_pos)) then + if (present(dgamma_dtau_pos)) then significantPositiveTau2: where(abs(tau_pos)-dst%threshold_stress(:,of) > tol_math_check) - dvel = 2.0_pReal*prm%burgers * prm%kink_height * prm%omega* effectiveLength & - * ( & - (needsGoodName + tau_pos * abs(needsGoodName)*BoltzmannRatio*prm%p & - * prm%q/prm%tau0 & - * StressRatio_pminus1*(1-StressRatio_p)**(prm%q-1.0_pReal) & - ) & - * ( 2.0_pReal*(prm%burgers**2.0_pReal)*tau_pos & - + prm%omega * prm%B* effectiveLength **2.0_pReal* needsGoodName & - ) & - - tau_pos * needsGoodName * (2.0_pReal*prm%burgers**2.0_pReal & - + prm%omega * prm%B *effectiveLength **2.0_pReal& - * (abs(needsGoodName)*BoltzmannRatio*prm%p *prm%q/prm%tau0 & - *StressRatio_pminus1*(1-StressRatio_p)**(prm%q-1.0_pReal) )& - ) & - ) & - /(2.0_pReal*prm%burgers**2.0_pReal*tau_pos & - + prm%omega * prm%B* effectiveLength**2.0_pReal* needsGoodName )**2.0_pReal + dtn = t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & + * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau0 + dtk = t_k / tau_pos + + dvel = prm%kink_height * (dtk + dtn) / (t_n + t_k)**2.0_pReal - dgdot_dtau_pos = DotGamma0 * dvel* 0.5_pReal + dgamma_dtau_pos = DotGamma0 * dvel* 0.5_pReal else where significantPositiveTau2 - dgdot_dtau_pos = 0.0_pReal + dgamma_dtau_pos = 0.0_pReal end where significantPositiveTau2 endif @@ -692,40 +680,27 @@ pure subroutine kinetics(Mp,Temperature,instance,of, & StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) - vel = 2.0_pReal*prm%burgers * prm%kink_height * prm%omega & - * effectiveLength * tau_neg * needsGoodName & - / ( 2.0_pReal*(prm%burgers**2.0_pReal)*tau_neg & - + prm%omega * prm%B * effectiveLength**2.0_pReal* needsGoodName & - ) + t_n = prm%burgers/(needsGoodName*prm%omega*effectiveLength) + t_k = effectiveLength * prm%B /(2.0_pReal*prm%burgers*tau_pos) ! our definition of tk is different with the one in dislotwin - gdot_neg = DotGamma0 * sign(vel,tau_neg) * 0.5_pReal + vel = prm%kink_height/(t_n + t_k) + + gamma_neg = DotGamma0 * sign(vel,tau_neg) * 0.5_pReal else where significantNegativeTau - gdot_neg = 0.0_pReal + gamma_neg = 0.0_pReal end where significantNegativeTau - if (present(dgdot_dtau_neg)) then + if (present(dgamma_dtau_neg)) then significantNegativeTau2: where(abs(tau_neg)-dst%threshold_stress(:,of) > tol_math_check) - dvel = 2.0_pReal*prm%burgers * prm%kink_height * prm%omega* effectiveLength & - * ( & - (needsGoodName + tau_neg * abs(needsGoodName)*BoltzmannRatio*prm%p & - * prm%q/prm%tau0 & - * StressRatio_pminus1*(1-StressRatio_p)**(prm%q-1.0_pReal) & - ) & - * ( 2.0_pReal*(prm%burgers**2.0_pReal)*tau_neg & - + prm%omega * prm%B* effectiveLength **2.0_pReal* needsGoodName & - ) & - - tau_neg * needsGoodName * (2.0_pReal*prm%burgers**2.0_pReal & - + prm%omega * prm%B *effectiveLength **2.0_pReal& - * (abs(needsGoodName)*BoltzmannRatio*prm%p *prm%q/prm%tau0 & - *StressRatio_pminus1*(1-StressRatio_p)**(prm%q-1.0_pReal) )& - ) & - ) & - /(2.0_pReal*prm%burgers**2.0_pReal*tau_neg & - + prm%omega * prm%B* effectiveLength**2.0_pReal* needsGoodName )**2.0_pReal + dtn = t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & + * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau0 + dtk = t_k / tau_pos + + dvel = prm%kink_height * (dtk + dtn) / (t_n + t_k)**2.0_pReal - dgdot_dtau_neg = DotGamma0 * dvel * 0.5_pReal + dgamma_dtau_neg = DotGamma0 * dvel * 0.5_pReal else where significantNegativeTau2 - dgdot_dtau_neg = 0.0_pReal + dgamma_dtau_neg = 0.0_pReal end where significantNegativeTau2 end if end associate From a12417a091d0b71e69e727685ea6b24a2282214b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 19 Mar 2019 05:46:47 +0100 Subject: [PATCH 04/12] as in dislotwin --- src/plastic_disloUCLA.f90 | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/src/plastic_disloUCLA.f90 b/src/plastic_disloUCLA.f90 index b099518ab..0d1d1814c 100644 --- a/src/plastic_disloUCLA.f90 +++ b/src/plastic_disloUCLA.f90 @@ -22,10 +22,10 @@ module plastic_disloUCLA enum, bind(c) enumerator :: & undefined_ID, & - rho_ID, & - rhoDip_ID, & - shearrate_ID, & - accumulatedshear_ID, & + rho_mob_ID, & + rho_dip_ID, & + gamma_dot_sl_ID, & + gamma_sl_ID, & mfp_ID, & thresholdstress_ID end enum @@ -292,13 +292,13 @@ subroutine plastic_disloUCLA_init() select case(trim(outputs(i))) case ('edge_density') - outputID = merge(rho_ID,undefined_ID,prm%totalNslip>0) + outputID = merge(rho_mob_ID,undefined_ID,prm%totalNslip>0) case ('dipole_density') - outputID = merge(rhoDip_ID,undefined_ID,prm%totalNslip>0) + outputID = merge(rho_dip_ID,undefined_ID,prm%totalNslip>0) case ('shear_rate','shearrate','shear_rate_slip','shearrate_slip') - outputID = merge(shearrate_ID,undefined_ID,prm%totalNslip>0) + outputID = merge(gamma_dot_sl_ID,undefined_ID,prm%totalNslip>0) case ('accumulated_shear','accumulatedshear','accumulated_shear_slip') - outputID = merge(accumulatedshear_ID,undefined_ID,prm%totalNslip>0) + outputID = merge(gamma_sl_ID,undefined_ID,prm%totalNslip>0) case ('mfp','mfp_slip') outputID = merge(mfp_ID,undefined_ID,prm%totalNslip>0) case ('threshold_stress','threshold_stress_slip') @@ -536,14 +536,14 @@ function plastic_disloUCLA_postResults(Mp,Temperature,instance,of) result(postRe outputsLoop: do o = 1,size(prm%outputID) select case(prm%outputID(o)) - case (rho_ID) + case (rho_mob_ID) postResults(c+1:c+prm%totalNslip) = stt%rho_mob(1:prm%totalNslip,of) - case (rhoDip_ID) + case (rho_dip_ID) postResults(c+1:c+prm%totalNslip) = stt%rho_dip(1:prm%totalNslip,of) - case (shearrate_ID) + case (gamma_dot_sl_ID) call kinetics(Mp,Temperature,instance,of,gdot_pos,gdot_neg) postResults(c+1:c+prm%totalNslip) = gdot_pos + gdot_neg - case (accumulatedshear_ID) + case (gamma_sl_ID) postResults(c+1:c+prm%totalNslip) = stt%gamma_sl(1:prm%totalNslip, of) case (mfp_ID) postResults(c+1:c+prm%totalNslip) = dst%mfp(1:prm%totalNslip, of) From 98d25ae48aaecc8c3723d5d9c4fb83e743733ced Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 19 Mar 2019 06:10:23 +0100 Subject: [PATCH 05/12] more name adjustments --- src/plastic_disloUCLA.f90 | 127 +++++++++++++++++++------------------- 1 file changed, 62 insertions(+), 65 deletions(-) diff --git a/src/plastic_disloUCLA.f90 b/src/plastic_disloUCLA.f90 index 0d1d1814c..27be6a022 100644 --- a/src/plastic_disloUCLA.f90 +++ b/src/plastic_disloUCLA.f90 @@ -26,7 +26,7 @@ module plastic_disloUCLA rho_dip_ID, & gamma_dot_sl_ID, & gamma_sl_ID, & - mfp_ID, & + Lambda_sl_ID, & thresholdstress_ID end enum @@ -34,22 +34,20 @@ module plastic_disloUCLA real(pReal) :: & aTolRho, & grainSize, & - SolidSolutionStrength, & !< Strength due to elements in solid solution mu, & D0, & !< prefactor for self-diffusion coefficient Qsd !< activation energy for dislocation climb real(pReal), dimension(:), allocatable :: & - rho_mob_0, & !< initial edge dislocation density + rho_mob_0, & !< initial edge dislocation density rhoDip0, & !< initial edge dipole density burgers, & !< absolute length of burgers vector [m] nonSchmidCoeff, & minDipDistance, & CLambda, & !< Adj. parameter for distance between 2 forest dislocations atomicVolume, & - tau_Peierls, & - tau0, & + tau_0, & !* mobility law parameters - H0kp, & !< activation energy for glide [J] + delta_F, & !< activation energy for glide [J] v0, & !< dislocation velocity prefactor [m/s] p, & !< p-exponent in glide velocity q, & !< q-exponent in glide velocity @@ -58,7 +56,7 @@ module plastic_disloUCLA w, & !< width of the kink pair omega !< attempt frequency for kink pair nucleation real(pReal), dimension(:,:), allocatable :: & - interaction_SlipSlip, & !< slip resistance from slip activity + h_sl_sl, & !< slip resistance from slip activity forestProjectionEdge real(pReal), dimension(:,:,:), allocatable :: & Schmid, & @@ -67,7 +65,7 @@ module plastic_disloUCLA integer :: & totalNslip !< total number of active slip system integer, dimension(:), allocatable :: & - Nslip !< number of active slip systems for each family + N_sl !< number of active slip systems for each family integer(kind(undefined_ID)), dimension(:),allocatable :: & outputID !< ID of each post result output logical :: & @@ -195,46 +193,45 @@ subroutine plastic_disloUCLA_init() !-------------------------------------------------------------------------------------------------- ! slip related parameters - prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) - prm%totalNslip = sum(prm%Nslip) + prm%N_sl = config%getInts('nslip',defaultVal=emptyIntArray) + prm%totalNslip = sum(prm%N_sl) slipActive: if (prm%totalNslip > 0) then - prm%Schmid = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),& + prm%Schmid = lattice_SchmidMatrix_slip(prm%N_sl,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) - prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1) + prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%N_sl,prm%nonSchmidCoeff,+1) + prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%N_sl,prm%nonSchmidCoeff,-1) else prm%nonSchmid_pos = prm%Schmid prm%nonSchmid_neg = prm%Schmid endif - prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(prm%Nslip, & - config%getFloats('interaction_slipslip'), & - config%getString('lattice_structure')) - prm%forestProjectionEdge = lattice_forestProjection(prm%Nslip,config%getString('lattice_structure'),& + prm%h_sl_sl = lattice_interaction_SlipBySlip(prm%N_sl, & + config%getFloats('interaction_slipslip'), & + config%getString('lattice_structure')) + prm%forestProjectionEdge = lattice_forestProjection(prm%N_sl,config%getString('lattice_structure'),& config%getFloat('c/a',defaultVal=0.0_pReal)) - prm%rho_mob_0 = config%getFloats('rhoedge0', requiredSize=size(prm%Nslip)) - prm%rhoDip0 = config%getFloats('rhoedgedip0', requiredSize=size(prm%Nslip)) - prm%v0 = config%getFloats('v0', requiredSize=size(prm%Nslip)) - prm%burgers = config%getFloats('slipburgers', requiredSize=size(prm%Nslip)) - prm%H0kp = config%getFloats('qedge', requiredSize=size(prm%Nslip)) + prm%rho_mob_0 = config%getFloats('rhoedge0', requiredSize=size(prm%N_sl)) + prm%rhoDip0 = config%getFloats('rhoedgedip0', requiredSize=size(prm%N_sl)) + prm%v0 = config%getFloats('v0', requiredSize=size(prm%N_sl)) + prm%burgers = config%getFloats('slipburgers', requiredSize=size(prm%N_sl)) + prm%delta_F = config%getFloats('qedge', requiredSize=size(prm%N_sl)) - prm%clambda = config%getFloats('clambdaslip', requiredSize=size(prm%Nslip)) - prm%tau_Peierls = config%getFloats('tau_peierls', requiredSize=size(prm%Nslip)) ! ToDo: Deprecated - prm%p = config%getFloats('p_slip', requiredSize=size(prm%Nslip), & - defaultVal=[(1.0_pReal,i=1,size(prm%Nslip))]) - prm%q = config%getFloats('q_slip', requiredSize=size(prm%Nslip), & - defaultVal=[(1.0_pReal,i=1,size(prm%Nslip))]) - prm%kink_height = config%getFloats('kink_height', requiredSize=size(prm%Nslip)) - prm%w = config%getFloats('kink_width', requiredSize=size(prm%Nslip)) - prm%omega = config%getFloats('omega', requiredSize=size(prm%Nslip)) - prm%B = config%getFloats('friction_coeff', requiredSize=size(prm%Nslip)) + prm%clambda = config%getFloats('clambdaslip', requiredSize=size(prm%N_sl)) + prm%tau_0 = config%getFloats('tau_peierls', requiredSize=size(prm%N_sl)) + prm%p = config%getFloats('p_slip', requiredSize=size(prm%N_sl), & + defaultVal=[(1.0_pReal,i=1,size(prm%N_sl))]) + prm%q = config%getFloats('q_slip', requiredSize=size(prm%N_sl), & + defaultVal=[(1.0_pReal,i=1,size(prm%N_sl))]) + prm%kink_height = config%getFloats('kink_height', requiredSize=size(prm%N_sl)) + prm%w = config%getFloats('kink_width', requiredSize=size(prm%N_sl)) + prm%omega = config%getFloats('omega', requiredSize=size(prm%N_sl)) + prm%B = config%getFloats('friction_coeff', requiredSize=size(prm%N_sl)) - prm%SolidSolutionStrength = config%getFloat('solidsolutionstrength') ! ToDo: Deprecated prm%grainSize = config%getFloat('grainsize') prm%D0 = config%getFloat('d0') prm%Qsd = config%getFloat('qsd') @@ -243,23 +240,22 @@ subroutine plastic_disloUCLA_init() prm%dipoleformation = config%getFloat('dipoleformationfactor') > 0.0_pReal !should be on by default, ToDo: change to /key/-type key ! expand: family => system - prm%rho_mob_0 = math_expand(prm%rho_mob_0, prm%Nslip) - prm%rhoDip0 = math_expand(prm%rhoDip0, prm%Nslip) - prm%q = math_expand(prm%q, prm%Nslip) - prm%p = math_expand(prm%p, prm%Nslip) - prm%H0kp = math_expand(prm%H0kp, prm%Nslip) - prm%burgers = math_expand(prm%burgers, prm%Nslip) - prm%kink_height = math_expand(prm%kink_height, prm%Nslip) - prm%w = math_expand(prm%w, prm%Nslip) - prm%omega = math_expand(prm%omega, prm%Nslip) - prm%tau_Peierls = math_expand(prm%tau_Peierls, prm%Nslip) - prm%v0 = math_expand(prm%v0, prm%Nslip) - prm%B = math_expand(prm%B, prm%Nslip) - prm%clambda = math_expand(prm%clambda, prm%Nslip) - prm%atomicVolume = math_expand(prm%atomicVolume, prm%Nslip) - prm%minDipDistance = math_expand(prm%minDipDistance, prm%Nslip) + prm%rho_mob_0 = math_expand(prm%rho_mob_0, prm%N_sl) + prm%rhoDip0 = math_expand(prm%rhoDip0, prm%N_sl) + prm%q = math_expand(prm%q, prm%N_sl) + prm%p = math_expand(prm%p, prm%N_sl) + prm%delta_F = math_expand(prm%delta_F, prm%N_sl) + prm%burgers = math_expand(prm%burgers, prm%N_sl) + prm%kink_height = math_expand(prm%kink_height, prm%N_sl) + prm%w = math_expand(prm%w, prm%N_sl) + prm%omega = math_expand(prm%omega, prm%N_sl) + prm%tau_0 = math_expand(prm%tau_0, prm%N_sl) + prm%v0 = math_expand(prm%v0, prm%N_sl) + prm%B = math_expand(prm%B, prm%N_sl) + prm%clambda = math_expand(prm%clambda, prm%N_sl) + prm%atomicVolume = math_expand(prm%atomicVolume, prm%N_sl) + prm%minDipDistance = math_expand(prm%minDipDistance, prm%N_sl) - prm%tau0 = prm%tau_peierls + prm%SolidSolutionStrength ! sanity checks if ( prm%D0 <= 0.0_pReal) extmsg = trim(extmsg)//' d0' @@ -268,8 +264,8 @@ subroutine plastic_disloUCLA_init() if (any(prm%rhoDip0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoedgedip0' if (any(prm%v0 < 0.0_pReal)) extmsg = trim(extmsg)//' v0' if (any(prm%burgers <= 0.0_pReal)) extmsg = trim(extmsg)//' slipburgers' - if (any(prm%H0kp <= 0.0_pReal)) extmsg = trim(extmsg)//' qedge' - if (any(prm%tau_peierls < 0.0_pReal)) extmsg = trim(extmsg)//' tau_peierls' + if (any(prm%delta_F <= 0.0_pReal)) extmsg = trim(extmsg)//' qedge' + if (any(prm%tau_0 < 0.0_pReal)) extmsg = trim(extmsg)//' tau_0' if (any(prm%minDipDistance <= 0.0_pReal)) extmsg = trim(extmsg)//' cedgedipmindistance or slipburgers' if (any(prm%atomicVolume <= 0.0_pReal)) extmsg = trim(extmsg)//' catomicvolume or slipburgers' @@ -300,7 +296,7 @@ subroutine plastic_disloUCLA_init() case ('accumulated_shear','accumulatedshear','accumulated_shear_slip') outputID = merge(gamma_sl_ID,undefined_ID,prm%totalNslip>0) case ('mfp','mfp_slip') - outputID = merge(mfp_ID,undefined_ID,prm%totalNslip>0) + outputID = merge(Lambda_sl_ID,undefined_ID,prm%totalNslip>0) case ('threshold_stress','threshold_stress_slip') outputID = merge(thresholdstress_ID,undefined_ID,prm%totalNslip>0) @@ -439,7 +435,7 @@ subroutine plastic_disloUCLA_dotState(Mp,Temperature,instance,of) call kinetics(Mp,Temperature,instance,of,& gdot_pos,gdot_neg, & - tau_pos1 = tau_pos,tau_neg1 = tau_neg) + tau_pos_out = tau_pos,tau_neg_out = tau_neg) dot%gamma_sl(:,of) = (gdot_pos+gdot_neg) ! ToDo: needs to be abs VacancyDiffusion = prm%D0*exp(-prm%Qsd/(kB*Temperature)) @@ -491,7 +487,7 @@ subroutine plastic_disloUCLA_dependentState(instance,of) prm%forestProjectionEdge(:,i))) dst%threshold_stress(i,of) = prm%mu*prm%burgers(i) & * sqrt(dot_product(stt%rho_mob(:,of)+stt%rho_dip(:,of), & - prm%interaction_SlipSlip(:,i))) + prm%h_sl_sl(:,i))) end forall dst%mfp(:,of) = prm%grainSize/(1.0_pReal+prm%grainSize*dst%dislocationSpacing(:,of)/prm%Clambda) @@ -545,7 +541,7 @@ function plastic_disloUCLA_postResults(Mp,Temperature,instance,of) result(postRe postResults(c+1:c+prm%totalNslip) = gdot_pos + gdot_neg case (gamma_sl_ID) postResults(c+1:c+prm%totalNslip) = stt%gamma_sl(1:prm%totalNslip, of) - case (mfp_ID) + case (Lambda_sl_ID) postResults(c+1:c+prm%totalNslip) = dst%mfp(1:prm%totalNslip, of) case (thresholdstress_ID) postResults(c+1:c+prm%totalNslip) = dst%threshold_stress(1:prm%totalNslip,of) @@ -595,7 +591,7 @@ end subroutine plastic_disloUCLA_results ! have the optional arguments at the end !-------------------------------------------------------------------------------------------------- pure subroutine kinetics(Mp,Temperature,instance,of, & - gamma_pos,gamma_neg,dgamma_dtau_pos,dgamma_dtau_neg,tau_pos1,tau_neg1) + gamma_pos,gamma_neg,dgamma_dtau_pos,dgamma_dtau_neg,tau_pos_out,tau_neg_out) use prec, only: & tol_math_check, & dEq, dNeq0 @@ -618,8 +614,8 @@ pure subroutine kinetics(Mp,Temperature,instance,of, & real(pReal), intent(out), optional, dimension(param(instance)%totalNslip) :: & dgamma_dtau_pos, & dgamma_dtau_neg, & - tau_pos1, & - tau_neg1 + tau_pos_out, & + tau_neg_out real(pReal), dimension(param(instance)%totalNslip) :: & StressRatio, & StressRatio_p,StressRatio_pminus1, & @@ -637,15 +633,15 @@ pure subroutine kinetics(Mp,Temperature,instance,of, & enddo - if (present(tau_pos1)) tau_pos1 = tau_pos - if (present(tau_neg1)) tau_neg1 = tau_neg + if (present(tau_pos_out)) tau_pos_out = tau_pos + if (present(tau_neg_out)) tau_neg_out = tau_neg - associate(BoltzmannRatio => prm%H0kp/(kB*Temperature), & + associate(BoltzmannRatio => prm%delta_F/(kB*Temperature), & DotGamma0 => stt%rho_mob(:,of)*prm%burgers*prm%v0, & effectiveLength => dst%mfp(:,of) - prm%w) significantPositiveTau: where(abs(tau_pos)-dst%threshold_stress(:,of) > tol_math_check) - StressRatio = (abs(tau_pos)-dst%threshold_stress(:,of))/prm%tau0 + StressRatio = (abs(tau_pos)-dst%threshold_stress(:,of))/prm%tau_0 StressRatio_p = StressRatio** prm%p StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) @@ -663,7 +659,7 @@ pure subroutine kinetics(Mp,Temperature,instance,of, & if (present(dgamma_dtau_pos)) then significantPositiveTau2: where(abs(tau_pos)-dst%threshold_stress(:,of) > tol_math_check) dtn = t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & - * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau0 + * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_0 dtk = t_k / tau_pos dvel = prm%kink_height * (dtk + dtn) / (t_n + t_k)**2.0_pReal @@ -675,7 +671,7 @@ pure subroutine kinetics(Mp,Temperature,instance,of, & endif significantNegativeTau: where(abs(tau_neg)-dst%threshold_stress(:,of) > tol_math_check) - StressRatio = (abs(tau_neg)-dst%threshold_stress(:,of))/prm%tau0 + StressRatio = (abs(tau_neg)-dst%threshold_stress(:,of))/prm%tau_0 StressRatio_p = StressRatio** prm%p StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) @@ -693,7 +689,7 @@ pure subroutine kinetics(Mp,Temperature,instance,of, & if (present(dgamma_dtau_neg)) then significantNegativeTau2: where(abs(tau_neg)-dst%threshold_stress(:,of) > tol_math_check) dtn = t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & - * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau0 + * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_0 dtk = t_k / tau_pos dvel = prm%kink_height * (dtk + dtn) / (t_n + t_k)**2.0_pReal @@ -703,6 +699,7 @@ pure subroutine kinetics(Mp,Temperature,instance,of, & dgamma_dtau_neg = 0.0_pReal end where significantNegativeTau2 end if + end associate end associate From 349acce9f28a860cec0b8b3d76fa80ed63f0016e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 20 Mar 2019 09:37:29 +0100 Subject: [PATCH 06/12] more renames to follow dislotwin and/or paper --- src/plastic_disloUCLA.f90 | 126 +++++++++++++++++++------------------- 1 file changed, 63 insertions(+), 63 deletions(-) diff --git a/src/plastic_disloUCLA.f90 b/src/plastic_disloUCLA.f90 index 27be6a022..2b9696863 100644 --- a/src/plastic_disloUCLA.f90 +++ b/src/plastic_disloUCLA.f90 @@ -33,14 +33,14 @@ module plastic_disloUCLA type, private :: tParameters real(pReal) :: & aTolRho, & - grainSize, & + D, & !< grain size mu, & D0, & !< prefactor for self-diffusion coefficient Qsd !< activation energy for dislocation climb real(pReal), dimension(:), allocatable :: & - rho_mob_0, & !< initial edge dislocation density - rhoDip0, & !< initial edge dipole density - burgers, & !< absolute length of burgers vector [m] + rho_mob_0, & !< initial dislocation density + rho_dip_0, & !< initial dipole density + b_sl, & !< magnitude of burgers vector [m] nonSchmidCoeff, & minDipDistance, & CLambda, & !< Adj. parameter for distance between 2 forest dislocations @@ -63,7 +63,7 @@ module plastic_disloUCLA nonSchmid_pos, & nonSchmid_neg integer :: & - totalNslip !< total number of active slip system + sum_N_sl !< total number of active slip system integer, dimension(:), allocatable :: & N_sl !< number of active slip systems for each family integer(kind(undefined_ID)), dimension(:),allocatable :: & @@ -194,8 +194,8 @@ subroutine plastic_disloUCLA_init() !-------------------------------------------------------------------------------------------------- ! slip related parameters prm%N_sl = config%getInts('nslip',defaultVal=emptyIntArray) - prm%totalNslip = sum(prm%N_sl) - slipActive: if (prm%totalNslip > 0) then + prm%sum_N_sl = sum(prm%N_sl) + slipActive: if (prm%sum_N_sl > 0) then prm%Schmid = lattice_SchmidMatrix_slip(prm%N_sl,config%getString('lattice_structure'),& config%getFloat('c/a',defaultVal=0.0_pReal)) @@ -216,9 +216,9 @@ subroutine plastic_disloUCLA_init() config%getFloat('c/a',defaultVal=0.0_pReal)) prm%rho_mob_0 = config%getFloats('rhoedge0', requiredSize=size(prm%N_sl)) - prm%rhoDip0 = config%getFloats('rhoedgedip0', requiredSize=size(prm%N_sl)) + prm%rho_dip_0 = config%getFloats('rhoedgedip0', requiredSize=size(prm%N_sl)) prm%v0 = config%getFloats('v0', requiredSize=size(prm%N_sl)) - prm%burgers = config%getFloats('slipburgers', requiredSize=size(prm%N_sl)) + prm%b_sl = config%getFloats('slipburgers', requiredSize=size(prm%N_sl)) prm%delta_F = config%getFloats('qedge', requiredSize=size(prm%N_sl)) prm%clambda = config%getFloats('clambdaslip', requiredSize=size(prm%N_sl)) @@ -232,20 +232,20 @@ subroutine plastic_disloUCLA_init() prm%omega = config%getFloats('omega', requiredSize=size(prm%N_sl)) prm%B = config%getFloats('friction_coeff', requiredSize=size(prm%N_sl)) - prm%grainSize = config%getFloat('grainsize') + prm%D = config%getFloat('grainsize') prm%D0 = config%getFloat('d0') prm%Qsd = config%getFloat('qsd') - prm%atomicVolume = config%getFloat('catomicvolume') * prm%burgers**3.0_pReal - prm%minDipDistance = config%getFloat('cedgedipmindistance') * prm%burgers + prm%atomicVolume = config%getFloat('catomicvolume') * prm%b_sl**3.0_pReal + prm%minDipDistance = config%getFloat('cedgedipmindistance') * prm%b_sl prm%dipoleformation = config%getFloat('dipoleformationfactor') > 0.0_pReal !should be on by default, ToDo: change to /key/-type key ! expand: family => system prm%rho_mob_0 = math_expand(prm%rho_mob_0, prm%N_sl) - prm%rhoDip0 = math_expand(prm%rhoDip0, prm%N_sl) + prm%rho_dip_0 = math_expand(prm%rho_dip_0, prm%N_sl) prm%q = math_expand(prm%q, prm%N_sl) prm%p = math_expand(prm%p, prm%N_sl) prm%delta_F = math_expand(prm%delta_F, prm%N_sl) - prm%burgers = math_expand(prm%burgers, prm%N_sl) + prm%b_sl = math_expand(prm%b_sl, prm%N_sl) prm%kink_height = math_expand(prm%kink_height, prm%N_sl) prm%w = math_expand(prm%w, prm%N_sl) prm%omega = math_expand(prm%omega, prm%N_sl) @@ -261,17 +261,17 @@ subroutine plastic_disloUCLA_init() if ( prm%D0 <= 0.0_pReal) extmsg = trim(extmsg)//' d0' if ( prm%Qsd <= 0.0_pReal) extmsg = trim(extmsg)//' qsd' if (any(prm%rho_mob_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoedge0' - if (any(prm%rhoDip0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoedgedip0' + if (any(prm%rho_dip_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoedgedip0' if (any(prm%v0 < 0.0_pReal)) extmsg = trim(extmsg)//' v0' - if (any(prm%burgers <= 0.0_pReal)) extmsg = trim(extmsg)//' slipburgers' + if (any(prm%b_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' slipb_sl' if (any(prm%delta_F <= 0.0_pReal)) extmsg = trim(extmsg)//' qedge' if (any(prm%tau_0 < 0.0_pReal)) extmsg = trim(extmsg)//' tau_0' - if (any(prm%minDipDistance <= 0.0_pReal)) extmsg = trim(extmsg)//' cedgedipmindistance or slipburgers' - if (any(prm%atomicVolume <= 0.0_pReal)) extmsg = trim(extmsg)//' catomicvolume or slipburgers' + if (any(prm%minDipDistance <= 0.0_pReal)) extmsg = trim(extmsg)//' cedgedipmindistance or slipb_sl' + if (any(prm%atomicVolume <= 0.0_pReal)) extmsg = trim(extmsg)//' catomicvolume or slipb_sl' else slipActive allocate(prm%rho_mob_0(0)) - allocate(prm%rhoDip0(0)) + allocate(prm%rho_dip_0(0)) endif slipActive !-------------------------------------------------------------------------------------------------- @@ -288,23 +288,23 @@ subroutine plastic_disloUCLA_init() select case(trim(outputs(i))) case ('edge_density') - outputID = merge(rho_mob_ID,undefined_ID,prm%totalNslip>0) + outputID = merge(rho_mob_ID,undefined_ID,prm%sum_N_sl>0) case ('dipole_density') - outputID = merge(rho_dip_ID,undefined_ID,prm%totalNslip>0) + outputID = merge(rho_dip_ID,undefined_ID,prm%sum_N_sl>0) case ('shear_rate','shearrate','shear_rate_slip','shearrate_slip') - outputID = merge(gamma_dot_sl_ID,undefined_ID,prm%totalNslip>0) + outputID = merge(gamma_dot_sl_ID,undefined_ID,prm%sum_N_sl>0) case ('accumulated_shear','accumulatedshear','accumulated_shear_slip') - outputID = merge(gamma_sl_ID,undefined_ID,prm%totalNslip>0) + outputID = merge(gamma_sl_ID,undefined_ID,prm%sum_N_sl>0) case ('mfp','mfp_slip') - outputID = merge(Lambda_sl_ID,undefined_ID,prm%totalNslip>0) + outputID = merge(Lambda_sl_ID,undefined_ID,prm%sum_N_sl>0) case ('threshold_stress','threshold_stress_slip') - outputID = merge(thresholdstress_ID,undefined_ID,prm%totalNslip>0) + outputID = merge(thresholdstress_ID,undefined_ID,prm%sum_N_sl>0) end select if (outputID /= undefined_ID) then plastic_disloUCLA_output(i,phase_plasticityInstance(p)) = outputs(i) - plastic_disloUCLA_sizePostResult(i,phase_plasticityInstance(p)) = prm%totalNslip + plastic_disloUCLA_sizePostResult(i,phase_plasticityInstance(p)) = prm%sum_N_sl prm%outputID = [prm%outputID, outputID] endif @@ -313,31 +313,31 @@ subroutine plastic_disloUCLA_init() !-------------------------------------------------------------------------------------------------- ! allocate state arrays NipcMyPhase = count(material_phase == p) - sizeDotState = size(['rhoEdge ','rhoEdgeDip ','accshearslip']) * prm%totalNslip + sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl sizeState = sizeDotState call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0, & - prm%totalNslip,0,0) + prm%sum_N_sl,0,0) plasticState(p)%sizePostResults = sum(plastic_disloUCLA_sizePostResult(:,phase_plasticityInstance(p))) !-------------------------------------------------------------------------------------------------- ! locally defined state aliases and initialization of state0 and aTolState startIndex = 1 - endIndex = prm%totalNslip + endIndex = prm%sum_N_sl stt%rho_mob=>plasticState(p)%state(startIndex:endIndex,:) stt%rho_mob= spread(prm%rho_mob_0,2,NipcMyPhase) dot%rho_mob=>plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolRho startIndex = endIndex + 1 - endIndex = endIndex + prm%totalNslip + endIndex = endIndex + prm%sum_N_sl stt%rho_dip=>plasticState(p)%state(startIndex:endIndex,:) - stt%rho_dip= spread(prm%rhoDip0,2,NipcMyPhase) + stt%rho_dip= spread(prm%rho_dip_0,2,NipcMyPhase) dot%rho_dip=>plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolRho startIndex = endIndex + 1 - endIndex = endIndex + prm%totalNslip + endIndex = endIndex + prm%sum_N_sl stt%gamma_sl=>plasticState(p)%state(startIndex:endIndex,:) dot%gamma_sl=>plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = 1.0e6_pReal !ToDo: better make optional parameter @@ -345,9 +345,9 @@ subroutine plastic_disloUCLA_init() plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%accumulatedSlip => plasticState(p)%state(startIndex:endIndex,:) - allocate(dst%mfp(prm%totalNslip,NipcMyPhase), source=0.0_pReal) - allocate(dst%dislocationSpacing(prm%totalNslip,NipcMyPhase),source=0.0_pReal) - allocate(dst%threshold_stress(prm%totalNslip,NipcMyPhase), source=0.0_pReal) + allocate(dst%mfp(prm%sum_N_sl,NipcMyPhase), source=0.0_pReal) + allocate(dst%dislocationSpacing(prm%sum_N_sl,NipcMyPhase),source=0.0_pReal) + allocate(dst%threshold_stress(prm%sum_N_sl,NipcMyPhase), source=0.0_pReal) plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally @@ -379,7 +379,7 @@ pure subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp,Mp,Temperature,inst integer :: & i,k,l,m,n - real(pReal), dimension(param(instance)%totalNslip) :: & + real(pReal), dimension(param(instance)%sum_N_sl) :: & gdot_pos,gdot_neg, & dgdot_dtau_pos,dgdot_dtau_neg @@ -389,7 +389,7 @@ pure subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp,Mp,Temperature,inst associate(prm => param(instance)) call kinetics(Mp,Temperature,instance,of,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) - do i = 1, prm%totalNslip + do i = 1, prm%sum_N_sl Lp = Lp + (gdot_pos(i)+gdot_neg(i))*prm%Schmid(1:3,1:3,i) forall (k=1:3,l=1:3,m=1:3,n=1:3) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & @@ -424,7 +424,7 @@ subroutine plastic_disloUCLA_dotState(Mp,Temperature,instance,of) real(pReal) :: & VacancyDiffusion - real(pReal), dimension(param(instance)%totalNslip) :: & + real(pReal), dimension(param(instance)%sum_N_sl) :: & gdot_pos, gdot_neg,& tau_pos,& tau_neg, & @@ -444,10 +444,10 @@ subroutine plastic_disloUCLA_dotState(Mp,Temperature,instance,of) DotRhoDipFormation = 0.0_pReal DotRhoEdgeDipClimb = 0.0_pReal else where - EdgeDipDistance = math_clip((3.0_pReal*prm%mu*prm%burgers)/(16.0_pReal*PI*abs(tau_pos)), & + EdgeDipDistance = math_clip((3.0_pReal*prm%mu*prm%b_sl)/(16.0_pReal*PI*abs(tau_pos)), & prm%minDipDistance, & ! lower limit dst%mfp(:,of)) ! upper limit - DotRhoDipFormation = merge(((2.0_pReal*EdgeDipDistance)/prm%burgers)* stt%rho_mob(:,of)*abs(dot%gamma_sl(:,of)), & ! ToDo: ignore region of spontaneous annihilation + DotRhoDipFormation = merge(((2.0_pReal*EdgeDipDistance)/prm%b_sl)* stt%rho_mob(:,of)*abs(dot%gamma_sl(:,of)), & ! ToDo: ignore region of spontaneous annihilation 0.0_pReal, & prm%dipoleformation) ClimbVelocity = (3.0_pReal*prm%mu*VacancyDiffusion*prm%atomicVolume/(2.0_pReal*pi*kB*Temperature)) & @@ -455,11 +455,11 @@ subroutine plastic_disloUCLA_dotState(Mp,Temperature,instance,of) DotRhoEdgeDipClimb = (4.0_pReal*ClimbVelocity*stt%rho_dip(:,of))/(EdgeDipDistance-prm%minDipDistance) ! ToDo: Discuss with Franz: Stress dependency? end where - dot%rho_mob(:,of) = abs(dot%gamma_sl(:,of))/(prm%burgers*dst%mfp(:,of)) & ! multiplication + dot%rho_mob(:,of) = abs(dot%gamma_sl(:,of))/(prm%b_sl*dst%mfp(:,of)) & ! multiplication - DotRhoDipFormation & - - (2.0_pReal*prm%minDipDistance)/prm%burgers*stt%rho_mob(:,of)*abs(dot%gamma_sl(:,of)) !* Spontaneous annihilation of 2 single edge dislocations + - (2.0_pReal*prm%minDipDistance)/prm%b_sl*stt%rho_mob(:,of)*abs(dot%gamma_sl(:,of)) !* Spontaneous annihilation of 2 single edge dislocations dot%rho_dip(:,of) = DotRhoDipFormation & - - (2.0_pReal*prm%minDipDistance)/prm%burgers*stt%rho_dip(:,of)*abs(dot%gamma_sl(:,of)) & !* Spontaneous annihilation of a single edge dislocation with a dipole constituent + - (2.0_pReal*prm%minDipDistance)/prm%b_sl*stt%rho_dip(:,of)*abs(dot%gamma_sl(:,of)) & !* Spontaneous annihilation of a single edge dislocation with a dipole constituent - DotRhoEdgeDipClimb end associate @@ -482,15 +482,15 @@ subroutine plastic_disloUCLA_dependentState(instance,of) associate(prm => param(instance), stt => state(instance),dst => dependentState(instance)) - forall (i = 1:prm%totalNslip) + forall (i = 1:prm%sum_N_sl) dst%dislocationSpacing(i,of) = sqrt(dot_product(stt%rho_mob(:,of)+stt%rho_dip(:,of), & prm%forestProjectionEdge(:,i))) - dst%threshold_stress(i,of) = prm%mu*prm%burgers(i) & + dst%threshold_stress(i,of) = prm%mu*prm%b_sl(i) & * sqrt(dot_product(stt%rho_mob(:,of)+stt%rho_dip(:,of), & prm%h_sl_sl(:,i))) end forall - dst%mfp(:,of) = prm%grainSize/(1.0_pReal+prm%grainSize*dst%dislocationSpacing(:,of)/prm%Clambda) + dst%mfp(:,of) = prm%D/(1.0_pReal+prm%D*dst%dislocationSpacing(:,of)/prm%Clambda) dst%dislocationSpacing(:,of) = dst%mfp(:,of) ! ToDo: Hack to recover wrong behavior for the moment end associate @@ -522,7 +522,7 @@ function plastic_disloUCLA_postResults(Mp,Temperature,instance,of) result(postRe integer :: & o,c,i - real(pReal), dimension(param(instance)%totalNslip) :: & + real(pReal), dimension(param(instance)%sum_N_sl) :: & gdot_pos,gdot_neg c = 0 @@ -533,22 +533,22 @@ function plastic_disloUCLA_postResults(Mp,Temperature,instance,of) result(postRe select case(prm%outputID(o)) case (rho_mob_ID) - postResults(c+1:c+prm%totalNslip) = stt%rho_mob(1:prm%totalNslip,of) + postResults(c+1:c+prm%sum_N_sl) = stt%rho_mob(1:prm%sum_N_sl,of) case (rho_dip_ID) - postResults(c+1:c+prm%totalNslip) = stt%rho_dip(1:prm%totalNslip,of) + postResults(c+1:c+prm%sum_N_sl) = stt%rho_dip(1:prm%sum_N_sl,of) case (gamma_dot_sl_ID) call kinetics(Mp,Temperature,instance,of,gdot_pos,gdot_neg) - postResults(c+1:c+prm%totalNslip) = gdot_pos + gdot_neg + postResults(c+1:c+prm%sum_N_sl) = gdot_pos + gdot_neg case (gamma_sl_ID) - postResults(c+1:c+prm%totalNslip) = stt%gamma_sl(1:prm%totalNslip, of) + postResults(c+1:c+prm%sum_N_sl) = stt%gamma_sl(1:prm%sum_N_sl, of) case (Lambda_sl_ID) - postResults(c+1:c+prm%totalNslip) = dst%mfp(1:prm%totalNslip, of) + postResults(c+1:c+prm%sum_N_sl) = dst%mfp(1:prm%sum_N_sl, of) case (thresholdstress_ID) - postResults(c+1:c+prm%totalNslip) = dst%threshold_stress(1:prm%totalNslip,of) + postResults(c+1:c+prm%sum_N_sl) = dst%threshold_stress(1:prm%sum_N_sl,of) end select - c = c + prm%totalNslip + c = c + prm%sum_N_sl enddo outputsLoop @@ -608,15 +608,15 @@ pure subroutine kinetics(Mp,Temperature,instance,of, & instance, & of - real(pReal), intent(out), dimension(param(instance)%totalNslip) :: & + real(pReal), intent(out), dimension(param(instance)%sum_N_sl) :: & gamma_pos, & gamma_neg - real(pReal), intent(out), optional, dimension(param(instance)%totalNslip) :: & + real(pReal), intent(out), optional, dimension(param(instance)%sum_N_sl) :: & dgamma_dtau_pos, & dgamma_dtau_neg, & tau_pos_out, & tau_neg_out - real(pReal), dimension(param(instance)%totalNslip) :: & + real(pReal), dimension(param(instance)%sum_N_sl) :: & StressRatio, & StressRatio_p,StressRatio_pminus1, & dvel, vel, & @@ -627,7 +627,7 @@ pure subroutine kinetics(Mp,Temperature,instance,of, & associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) - do j = 1, prm%totalNslip + do j = 1, prm%sum_N_sl tau_pos(j) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,j)) tau_neg(j) = math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,j)) enddo @@ -637,7 +637,7 @@ pure subroutine kinetics(Mp,Temperature,instance,of, & if (present(tau_neg_out)) tau_neg_out = tau_neg associate(BoltzmannRatio => prm%delta_F/(kB*Temperature), & - DotGamma0 => stt%rho_mob(:,of)*prm%burgers*prm%v0, & + DotGamma0 => stt%rho_mob(:,of)*prm%b_sl*prm%v0, & effectiveLength => dst%mfp(:,of) - prm%w) significantPositiveTau: where(abs(tau_pos)-dst%threshold_stress(:,of) > tol_math_check) @@ -646,8 +646,8 @@ pure subroutine kinetics(Mp,Temperature,instance,of, & StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) - t_n = prm%burgers/(needsGoodName*prm%omega*effectiveLength) - t_k = effectiveLength * prm%B /(2.0_pReal*prm%burgers*tau_pos) ! our definition of tk is different with the one in dislotwin + t_n = prm%b_sl/(needsGoodName*prm%omega*effectiveLength) + t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos) ! our definition of tk is different with the one in dislotwin vel = prm%kink_height/(t_n + t_k) @@ -676,8 +676,8 @@ pure subroutine kinetics(Mp,Temperature,instance,of, & StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) - t_n = prm%burgers/(needsGoodName*prm%omega*effectiveLength) - t_k = effectiveLength * prm%B /(2.0_pReal*prm%burgers*tau_pos) ! our definition of tk is different with the one in dislotwin + t_n = prm%b_sl/(needsGoodName*prm%omega*effectiveLength) + t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos) ! our definition of tk is different with the one in dislotwin vel = prm%kink_height/(t_n + t_k) From 97f6566d989b33456de41f8f728ce93b349dc951 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 20 Mar 2019 16:52:43 +0100 Subject: [PATCH 07/12] more renames --- src/plastic_disloUCLA.f90 | 75 +++++++++++++++++++-------------------- 1 file changed, 37 insertions(+), 38 deletions(-) diff --git a/src/plastic_disloUCLA.f90 b/src/plastic_disloUCLA.f90 index 2b9696863..b9fc29ad1 100644 --- a/src/plastic_disloUCLA.f90 +++ b/src/plastic_disloUCLA.f90 @@ -36,14 +36,14 @@ module plastic_disloUCLA D, & !< grain size mu, & D0, & !< prefactor for self-diffusion coefficient - Qsd !< activation energy for dislocation climb + Q_cl !< activation energy for dislocation climb real(pReal), dimension(:), allocatable :: & rho_mob_0, & !< initial dislocation density rho_dip_0, & !< initial dipole density b_sl, & !< magnitude of burgers vector [m] nonSchmidCoeff, & - minDipDistance, & - CLambda, & !< Adj. parameter for distance between 2 forest dislocations + D_a, & + i_sl, & !< Adj. parameter for distance between 2 forest dislocations atomicVolume, & tau_0, & !* mobility law parameters @@ -81,8 +81,7 @@ module plastic_disloUCLA type, private :: tDisloUCLAdependentState real(pReal), dimension(:,:), allocatable :: & - mfp, & - dislocationSpacing, & + Lambda_sl, & threshold_stress end type tDisloUCLAdependentState @@ -221,7 +220,7 @@ subroutine plastic_disloUCLA_init() prm%b_sl = config%getFloats('slipburgers', requiredSize=size(prm%N_sl)) prm%delta_F = config%getFloats('qedge', requiredSize=size(prm%N_sl)) - prm%clambda = config%getFloats('clambdaslip', requiredSize=size(prm%N_sl)) + prm%i_sl = config%getFloats('clambdaslip', requiredSize=size(prm%N_sl)) prm%tau_0 = config%getFloats('tau_peierls', requiredSize=size(prm%N_sl)) prm%p = config%getFloats('p_slip', requiredSize=size(prm%N_sl), & defaultVal=[(1.0_pReal,i=1,size(prm%N_sl))]) @@ -232,12 +231,12 @@ subroutine plastic_disloUCLA_init() prm%omega = config%getFloats('omega', requiredSize=size(prm%N_sl)) prm%B = config%getFloats('friction_coeff', requiredSize=size(prm%N_sl)) - prm%D = config%getFloat('grainsize') - prm%D0 = config%getFloat('d0') - prm%Qsd = config%getFloat('qsd') - prm%atomicVolume = config%getFloat('catomicvolume') * prm%b_sl**3.0_pReal - prm%minDipDistance = config%getFloat('cedgedipmindistance') * prm%b_sl - prm%dipoleformation = config%getFloat('dipoleformationfactor') > 0.0_pReal !should be on by default, ToDo: change to /key/-type key + prm%D = config%getFloat('grainsize') + prm%D0 = config%getFloat('d0') + prm%Q_cl = config%getFloat('qsd') + prm%atomicVolume = config%getFloat('catomicvolume') * prm%b_sl**3.0_pReal + prm%D_a = config%getFloat('cedgedipmindistance') * prm%b_sl + prm%dipoleformation = config%getFloat('dipoleformationfactor') > 0.0_pReal !should be on by default, ToDo: change to /key/-type key ! expand: family => system prm%rho_mob_0 = math_expand(prm%rho_mob_0, prm%N_sl) @@ -252,21 +251,21 @@ subroutine plastic_disloUCLA_init() prm%tau_0 = math_expand(prm%tau_0, prm%N_sl) prm%v0 = math_expand(prm%v0, prm%N_sl) prm%B = math_expand(prm%B, prm%N_sl) - prm%clambda = math_expand(prm%clambda, prm%N_sl) + prm%i_sl = math_expand(prm%i_sl, prm%N_sl) prm%atomicVolume = math_expand(prm%atomicVolume, prm%N_sl) - prm%minDipDistance = math_expand(prm%minDipDistance, prm%N_sl) + prm%D_a = math_expand(prm%D_a, prm%N_sl) ! sanity checks if ( prm%D0 <= 0.0_pReal) extmsg = trim(extmsg)//' d0' - if ( prm%Qsd <= 0.0_pReal) extmsg = trim(extmsg)//' qsd' + if ( prm%Q_cl <= 0.0_pReal) extmsg = trim(extmsg)//' Q_cl' if (any(prm%rho_mob_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoedge0' if (any(prm%rho_dip_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoedgedip0' if (any(prm%v0 < 0.0_pReal)) extmsg = trim(extmsg)//' v0' if (any(prm%b_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' slipb_sl' if (any(prm%delta_F <= 0.0_pReal)) extmsg = trim(extmsg)//' qedge' if (any(prm%tau_0 < 0.0_pReal)) extmsg = trim(extmsg)//' tau_0' - if (any(prm%minDipDistance <= 0.0_pReal)) extmsg = trim(extmsg)//' cedgedipmindistance or slipb_sl' + if (any(prm%D_a <= 0.0_pReal)) extmsg = trim(extmsg)//' cedgedipmindistance or slipb_sl' if (any(prm%atomicVolume <= 0.0_pReal)) extmsg = trim(extmsg)//' catomicvolume or slipb_sl' else slipActive @@ -345,8 +344,7 @@ subroutine plastic_disloUCLA_init() plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%accumulatedSlip => plasticState(p)%state(startIndex:endIndex,:) - allocate(dst%mfp(prm%sum_N_sl,NipcMyPhase), source=0.0_pReal) - allocate(dst%dislocationSpacing(prm%sum_N_sl,NipcMyPhase),source=0.0_pReal) + allocate(dst%Lambda_sl(prm%sum_N_sl,NipcMyPhase), source=0.0_pReal) allocate(dst%threshold_stress(prm%sum_N_sl,NipcMyPhase), source=0.0_pReal) plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally @@ -405,7 +403,7 @@ end subroutine plastic_disloUCLA_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -subroutine plastic_disloUCLA_dotState(Mp,Temperature,instance,of) +subroutine plastic_disloUCLA_dotState(Mp,T,instance,of) use prec, only: & tol_math_check, & dEq0 @@ -417,7 +415,7 @@ subroutine plastic_disloUCLA_dotState(Mp,Temperature,instance,of) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress real(pReal), intent(in) :: & - temperature !< temperature + T !< temperature integer, intent(in) :: & instance, & of @@ -428,39 +426,39 @@ subroutine plastic_disloUCLA_dotState(Mp,Temperature,instance,of) gdot_pos, gdot_neg,& tau_pos,& tau_neg, & - DotRhoDipFormation, ClimbVelocity, EdgeDipDistance, & + DotRhoDipFormation, v_cl, EdgeDipDistance, & DotRhoEdgeDipClimb associate(prm => param(instance), stt => state(instance),dot => dotState(instance), dst => dependentState(instance)) - call kinetics(Mp,Temperature,instance,of,& + call kinetics(Mp,T,instance,of,& gdot_pos,gdot_neg, & tau_pos_out = tau_pos,tau_neg_out = tau_neg) dot%gamma_sl(:,of) = (gdot_pos+gdot_neg) ! ToDo: needs to be abs - VacancyDiffusion = prm%D0*exp(-prm%Qsd/(kB*Temperature)) + VacancyDiffusion = prm%D0*exp(-prm%Q_cl/(kB*T)) where(dEq0(tau_pos)) ! ToDo: use avg of pos and neg DotRhoDipFormation = 0.0_pReal DotRhoEdgeDipClimb = 0.0_pReal else where EdgeDipDistance = math_clip((3.0_pReal*prm%mu*prm%b_sl)/(16.0_pReal*PI*abs(tau_pos)), & - prm%minDipDistance, & ! lower limit - dst%mfp(:,of)) ! upper limit - DotRhoDipFormation = merge(((2.0_pReal*EdgeDipDistance)/prm%b_sl)* stt%rho_mob(:,of)*abs(dot%gamma_sl(:,of)), & ! ToDo: ignore region of spontaneous annihilation + prm%D_a, & ! lower limit + dst%Lambda_sl(:,of)) ! upper limit + dotRhoDipFormation = merge(((2.0_pReal*EdgeDipDistance)/prm%b_sl)* stt%rho_mob(:,of)*abs(dot%gamma_sl(:,of)), & ! ToDo: ignore region of spontaneous annihilation 0.0_pReal, & prm%dipoleformation) - ClimbVelocity = (3.0_pReal*prm%mu*VacancyDiffusion*prm%atomicVolume/(2.0_pReal*pi*kB*Temperature)) & - * (1.0_pReal/(EdgeDipDistance+prm%minDipDistance)) - DotRhoEdgeDipClimb = (4.0_pReal*ClimbVelocity*stt%rho_dip(:,of))/(EdgeDipDistance-prm%minDipDistance) ! ToDo: Discuss with Franz: Stress dependency? + v_cl = (3.0_pReal*prm%mu*VacancyDiffusion*prm%atomicVolume/(2.0_pReal*pi*kB*T)) & + * (1.0_pReal/(EdgeDipDistance+prm%D_a)) + dotRhoEdgeDipClimb = (4.0_pReal*v_cl*stt%rho_dip(:,of))/(EdgeDipDistance-prm%D_a) ! ToDo: Discuss with Franz: Stress dependency? end where - dot%rho_mob(:,of) = abs(dot%gamma_sl(:,of))/(prm%b_sl*dst%mfp(:,of)) & ! multiplication + dot%rho_mob(:,of) = abs(dot%gamma_sl(:,of))/(prm%b_sl*dst%Lambda_sl(:,of)) & ! multiplication - DotRhoDipFormation & - - (2.0_pReal*prm%minDipDistance)/prm%b_sl*stt%rho_mob(:,of)*abs(dot%gamma_sl(:,of)) !* Spontaneous annihilation of 2 single edge dislocations + - (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_mob(:,of)*abs(dot%gamma_sl(:,of)) ! Spontaneous annihilation of 2 single edge dislocations dot%rho_dip(:,of) = DotRhoDipFormation & - - (2.0_pReal*prm%minDipDistance)/prm%b_sl*stt%rho_dip(:,of)*abs(dot%gamma_sl(:,of)) & !* Spontaneous annihilation of a single edge dislocation with a dipole constituent - - DotRhoEdgeDipClimb + - (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_dip(:,of)*abs(dot%gamma_sl(:,of)) & ! Spontaneous annihilation of a single edge dislocation with a dipole constituent + - DotRhoEdgeDipClimb end associate @@ -477,21 +475,22 @@ subroutine plastic_disloUCLA_dependentState(instance,of) instance, & of + real(pReal), dimension(param(instance)%sum_N_sl) :: & + dislocationSpacing integer :: & i associate(prm => param(instance), stt => state(instance),dst => dependentState(instance)) forall (i = 1:prm%sum_N_sl) - dst%dislocationSpacing(i,of) = sqrt(dot_product(stt%rho_mob(:,of)+stt%rho_dip(:,of), & + dislocationSpacing(i) = sqrt(dot_product(stt%rho_mob(:,of)+stt%rho_dip(:,of), & prm%forestProjectionEdge(:,i))) dst%threshold_stress(i,of) = prm%mu*prm%b_sl(i) & * sqrt(dot_product(stt%rho_mob(:,of)+stt%rho_dip(:,of), & prm%h_sl_sl(:,i))) end forall - dst%mfp(:,of) = prm%D/(1.0_pReal+prm%D*dst%dislocationSpacing(:,of)/prm%Clambda) - dst%dislocationSpacing(:,of) = dst%mfp(:,of) ! ToDo: Hack to recover wrong behavior for the moment + dst%Lambda_sl(:,of) = prm%D/(1.0_pReal+prm%D*dislocationSpacing/prm%i_sl) end associate @@ -542,7 +541,7 @@ function plastic_disloUCLA_postResults(Mp,Temperature,instance,of) result(postRe case (gamma_sl_ID) postResults(c+1:c+prm%sum_N_sl) = stt%gamma_sl(1:prm%sum_N_sl, of) case (Lambda_sl_ID) - postResults(c+1:c+prm%sum_N_sl) = dst%mfp(1:prm%sum_N_sl, of) + postResults(c+1:c+prm%sum_N_sl) = dst%Lambda_sl(1:prm%sum_N_sl, of) case (thresholdstress_ID) postResults(c+1:c+prm%sum_N_sl) = dst%threshold_stress(1:prm%sum_N_sl,of) @@ -638,7 +637,7 @@ pure subroutine kinetics(Mp,Temperature,instance,of, & associate(BoltzmannRatio => prm%delta_F/(kB*Temperature), & DotGamma0 => stt%rho_mob(:,of)*prm%b_sl*prm%v0, & - effectiveLength => dst%mfp(:,of) - prm%w) + effectiveLength => dst%Lambda_sl(:,of) - prm%w) significantPositiveTau: where(abs(tau_pos)-dst%threshold_stress(:,of) > tol_math_check) StressRatio = (abs(tau_pos)-dst%threshold_stress(:,of))/prm%tau_0 From 25f2d78656347e94ae87254ecb57fb22e8539801 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 21 Mar 2019 00:06:31 +0100 Subject: [PATCH 08/12] final renames (for the moment) --- src/plastic_disloUCLA.f90 | 48 +++++++++++++++++++-------------------- 1 file changed, 24 insertions(+), 24 deletions(-) diff --git a/src/plastic_disloUCLA.f90 b/src/plastic_disloUCLA.f90 index b9fc29ad1..9dfc63c01 100644 --- a/src/plastic_disloUCLA.f90 +++ b/src/plastic_disloUCLA.f90 @@ -24,7 +24,7 @@ module plastic_disloUCLA undefined_ID, & rho_mob_ID, & rho_dip_ID, & - gamma_dot_sl_ID, & + dot_gamma_sl_ID, & gamma_sl_ID, & Lambda_sl_ID, & thresholdstress_ID @@ -32,10 +32,10 @@ module plastic_disloUCLA type, private :: tParameters real(pReal) :: & - aTolRho, & + aTol_rho, & D, & !< grain size mu, & - D0, & !< prefactor for self-diffusion coefficient + D_0, & !< prefactor for self-diffusion coefficient Q_cl !< activation energy for dislocation climb real(pReal), dimension(:), allocatable :: & rho_mob_0, & !< initial dislocation density @@ -185,10 +185,10 @@ subroutine plastic_disloUCLA_init() ! optional parameters that need to be defined prm%mu = lattice_mu(p) - prm%aTolRho = config%getFloat('atol_rho') + prm%aTol_rho = config%getFloat('atol_rho') ! sanity checks - if (prm%aTolRho <= 0.0_pReal) extmsg = trim(extmsg)//' atol_rho' + if (prm%aTol_rho <= 0.0_pReal) extmsg = trim(extmsg)//' atol_rho' !-------------------------------------------------------------------------------------------------- ! slip related parameters @@ -232,7 +232,7 @@ subroutine plastic_disloUCLA_init() prm%B = config%getFloats('friction_coeff', requiredSize=size(prm%N_sl)) prm%D = config%getFloat('grainsize') - prm%D0 = config%getFloat('d0') + prm%D_0 = config%getFloat('d0') prm%Q_cl = config%getFloat('qsd') prm%atomicVolume = config%getFloat('catomicvolume') * prm%b_sl**3.0_pReal prm%D_a = config%getFloat('cedgedipmindistance') * prm%b_sl @@ -257,7 +257,7 @@ subroutine plastic_disloUCLA_init() ! sanity checks - if ( prm%D0 <= 0.0_pReal) extmsg = trim(extmsg)//' d0' + if ( prm%D_0 <= 0.0_pReal) extmsg = trim(extmsg)//' D_0' if ( prm%Q_cl <= 0.0_pReal) extmsg = trim(extmsg)//' Q_cl' if (any(prm%rho_mob_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoedge0' if (any(prm%rho_dip_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoedgedip0' @@ -291,7 +291,7 @@ subroutine plastic_disloUCLA_init() case ('dipole_density') outputID = merge(rho_dip_ID,undefined_ID,prm%sum_N_sl>0) case ('shear_rate','shearrate','shear_rate_slip','shearrate_slip') - outputID = merge(gamma_dot_sl_ID,undefined_ID,prm%sum_N_sl>0) + outputID = merge(dot_gamma_sl_ID,undefined_ID,prm%sum_N_sl>0) case ('accumulated_shear','accumulatedshear','accumulated_shear_slip') outputID = merge(gamma_sl_ID,undefined_ID,prm%sum_N_sl>0) case ('mfp','mfp_slip') @@ -326,20 +326,20 @@ subroutine plastic_disloUCLA_init() stt%rho_mob=>plasticState(p)%state(startIndex:endIndex,:) stt%rho_mob= spread(prm%rho_mob_0,2,NipcMyPhase) dot%rho_mob=>plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolRho + plasticState(p)%aTolState(startIndex:endIndex) = prm%aTol_rho startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl stt%rho_dip=>plasticState(p)%state(startIndex:endIndex,:) stt%rho_dip= spread(prm%rho_dip_0,2,NipcMyPhase) dot%rho_dip=>plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolRho + plasticState(p)%aTolState(startIndex:endIndex) = prm%aTol_rho startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl stt%gamma_sl=>plasticState(p)%state(startIndex:endIndex,:) dot%gamma_sl=>plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%aTolState(startIndex:endIndex) = 1.0e6_pReal !ToDo: better make optional parameter + plasticState(p)%aTolState(startIndex:endIndex) = 1.0e6_pReal !Don't use for convergence check ! global alias plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%accumulatedSlip => plasticState(p)%state(startIndex:endIndex,:) @@ -359,7 +359,7 @@ end subroutine plastic_disloUCLA_init !-------------------------------------------------------------------------------------------------- !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- -pure subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp,Mp,Temperature,instance,of) +pure subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) implicit none real(pReal), dimension(3,3), intent(out) :: & @@ -370,7 +370,7 @@ pure subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp,Mp,Temperature,inst real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress real(pReal), intent(in) :: & - temperature !< temperature + T !< temperature integer, intent(in) :: & instance, & of @@ -386,7 +386,7 @@ pure subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp,Mp,Temperature,inst associate(prm => param(instance)) - call kinetics(Mp,Temperature,instance,of,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) + call kinetics(Mp,T,instance,of,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) do i = 1, prm%sum_N_sl Lp = Lp + (gdot_pos(i)+gdot_neg(i))*prm%Schmid(1:3,1:3,i) forall (k=1:3,l=1:3,m=1:3,n=1:3) & @@ -435,10 +435,10 @@ subroutine plastic_disloUCLA_dotState(Mp,T,instance,of) gdot_pos,gdot_neg, & tau_pos_out = tau_pos,tau_neg_out = tau_neg) - dot%gamma_sl(:,of) = (gdot_pos+gdot_neg) ! ToDo: needs to be abs - VacancyDiffusion = prm%D0*exp(-prm%Q_cl/(kB*T)) + dot%gamma_sl(:,of) = (gdot_pos+gdot_neg) ! ToDo: needs to be abs + VacancyDiffusion = prm%D_0*exp(-prm%Q_cl/(kB*T)) - where(dEq0(tau_pos)) ! ToDo: use avg of pos and neg + where(dEq0(tau_pos)) ! ToDo: use avg of pos and neg DotRhoDipFormation = 0.0_pReal DotRhoEdgeDipClimb = 0.0_pReal else where @@ -500,7 +500,7 @@ end subroutine plastic_disloUCLA_dependentState !-------------------------------------------------------------------------------------------------- !> @brief return array of constitutive results !-------------------------------------------------------------------------------------------------- -function plastic_disloUCLA_postResults(Mp,Temperature,instance,of) result(postResults) +function plastic_disloUCLA_postResults(Mp,T,instance,of) result(postResults) use prec, only: & dEq, dNeq0 use math, only: & @@ -511,7 +511,7 @@ function plastic_disloUCLA_postResults(Mp,Temperature,instance,of) result(postRe real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress real(pReal), intent(in) :: & - temperature !< temperature + T !< temperature integer, intent(in) :: & instance, & of @@ -535,8 +535,8 @@ function plastic_disloUCLA_postResults(Mp,Temperature,instance,of) result(postRe postResults(c+1:c+prm%sum_N_sl) = stt%rho_mob(1:prm%sum_N_sl,of) case (rho_dip_ID) postResults(c+1:c+prm%sum_N_sl) = stt%rho_dip(1:prm%sum_N_sl,of) - case (gamma_dot_sl_ID) - call kinetics(Mp,Temperature,instance,of,gdot_pos,gdot_neg) + case (dot_gamma_sl_ID) + call kinetics(Mp,T,instance,of,gdot_pos,gdot_neg) postResults(c+1:c+prm%sum_N_sl) = gdot_pos + gdot_neg case (gamma_sl_ID) postResults(c+1:c+prm%sum_N_sl) = stt%gamma_sl(1:prm%sum_N_sl, of) @@ -589,7 +589,7 @@ end subroutine plastic_disloUCLA_results ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end !-------------------------------------------------------------------------------------------------- -pure subroutine kinetics(Mp,Temperature,instance,of, & +pure subroutine kinetics(Mp,T,instance,of, & gamma_pos,gamma_neg,dgamma_dtau_pos,dgamma_dtau_neg,tau_pos_out,tau_neg_out) use prec, only: & tol_math_check, & @@ -602,7 +602,7 @@ pure subroutine kinetics(Mp,Temperature,instance,of, & real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress real(pReal), intent(in) :: & - temperature !< temperature + T !< temperature integer, intent(in) :: & instance, & of @@ -635,7 +635,7 @@ pure subroutine kinetics(Mp,Temperature,instance,of, & if (present(tau_pos_out)) tau_pos_out = tau_pos if (present(tau_neg_out)) tau_neg_out = tau_neg - associate(BoltzmannRatio => prm%delta_F/(kB*Temperature), & + associate(BoltzmannRatio => prm%delta_F/(kB*T), & DotGamma0 => stt%rho_mob(:,of)*prm%b_sl*prm%v0, & effectiveLength => dst%Lambda_sl(:,of) - prm%w) From fa9604232065dfcedf712236b31cac555c71a8bb Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 22 Mar 2019 10:49:24 +0100 Subject: [PATCH 09/12] correct names for shear rates --- src/plastic_disloUCLA.f90 | 46 +++++++++++++++++++-------------------- 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/src/plastic_disloUCLA.f90 b/src/plastic_disloUCLA.f90 index 9dfc63c01..22bb8ff96 100644 --- a/src/plastic_disloUCLA.f90 +++ b/src/plastic_disloUCLA.f90 @@ -378,21 +378,21 @@ pure subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) integer :: & i,k,l,m,n real(pReal), dimension(param(instance)%sum_N_sl) :: & - gdot_pos,gdot_neg, & - dgdot_dtau_pos,dgdot_dtau_neg + dot_gamma_pos,dot_gamma_neg, & + ddot_gamma_dtau_pos,ddot_gamma_dtau_neg Lp = 0.0_pReal dLp_dMp = 0.0_pReal associate(prm => param(instance)) - call kinetics(Mp,T,instance,of,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) + call kinetics(Mp,T,instance,of,dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg) do i = 1, prm%sum_N_sl - Lp = Lp + (gdot_pos(i)+gdot_neg(i))*prm%Schmid(1:3,1:3,i) + Lp = Lp + (dot_gamma_pos(i)+dot_gamma_neg(i))*prm%Schmid(1:3,1:3,i) forall (k=1:3,l=1:3,m=1:3,n=1:3) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & - + dgdot_dtau_pos(i) * prm%Schmid(k,l,i) * prm%nonSchmid_pos(m,n,i) & - + dgdot_dtau_neg(i) * prm%Schmid(k,l,i) * prm%nonSchmid_neg(m,n,i) + + ddot_gamma_dtau_pos(i) * prm%Schmid(k,l,i) * prm%nonSchmid_pos(m,n,i) & + + ddot_gamma_dtau_neg(i) * prm%Schmid(k,l,i) * prm%nonSchmid_neg(m,n,i) enddo end associate @@ -590,7 +590,7 @@ end subroutine plastic_disloUCLA_results ! have the optional arguments at the end !-------------------------------------------------------------------------------------------------- pure subroutine kinetics(Mp,T,instance,of, & - gamma_pos,gamma_neg,dgamma_dtau_pos,dgamma_dtau_neg,tau_pos_out,tau_neg_out) + dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg,tau_pos_out,tau_neg_out) use prec, only: & tol_math_check, & dEq, dNeq0 @@ -608,11 +608,11 @@ pure subroutine kinetics(Mp,T,instance,of, & of real(pReal), intent(out), dimension(param(instance)%sum_N_sl) :: & - gamma_pos, & - gamma_neg + dot_gamma_pos, & + dot_gamma_neg real(pReal), intent(out), optional, dimension(param(instance)%sum_N_sl) :: & - dgamma_dtau_pos, & - dgamma_dtau_neg, & + ddot_gamma_dtau_pos, & + ddot_gamma_dtau_neg, & tau_pos_out, & tau_neg_out real(pReal), dimension(param(instance)%sum_N_sl) :: & @@ -636,7 +636,7 @@ pure subroutine kinetics(Mp,T,instance,of, & if (present(tau_neg_out)) tau_neg_out = tau_neg associate(BoltzmannRatio => prm%delta_F/(kB*T), & - DotGamma0 => stt%rho_mob(:,of)*prm%b_sl*prm%v0, & + dot_gamma_0 => stt%rho_mob(:,of)*prm%b_sl*prm%v0, & effectiveLength => dst%Lambda_sl(:,of) - prm%w) significantPositiveTau: where(abs(tau_pos)-dst%threshold_stress(:,of) > tol_math_check) @@ -650,12 +650,12 @@ pure subroutine kinetics(Mp,T,instance,of, & vel = prm%kink_height/(t_n + t_k) - gamma_pos = DotGamma0 * sign(vel,tau_pos) * 0.5_pReal + dot_gamma_pos = dot_gamma_0 * sign(vel,tau_pos) * 0.5_pReal else where significantPositiveTau - gamma_pos = 0.0_pReal + dot_gamma_pos = 0.0_pReal end where significantPositiveTau - if (present(dgamma_dtau_pos)) then + if (present(ddot_gamma_dtau_pos)) then significantPositiveTau2: where(abs(tau_pos)-dst%threshold_stress(:,of) > tol_math_check) dtn = t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_0 @@ -663,9 +663,9 @@ pure subroutine kinetics(Mp,T,instance,of, & dvel = prm%kink_height * (dtk + dtn) / (t_n + t_k)**2.0_pReal - dgamma_dtau_pos = DotGamma0 * dvel* 0.5_pReal + ddot_gamma_dtau_pos = dot_gamma_0 * dvel* 0.5_pReal else where significantPositiveTau2 - dgamma_dtau_pos = 0.0_pReal + ddot_gamma_dtau_pos = 0.0_pReal end where significantPositiveTau2 endif @@ -680,22 +680,22 @@ pure subroutine kinetics(Mp,T,instance,of, & vel = prm%kink_height/(t_n + t_k) - gamma_neg = DotGamma0 * sign(vel,tau_neg) * 0.5_pReal + dot_gamma_neg = dot_gamma_0 * sign(vel,tau_neg) * 0.5_pReal else where significantNegativeTau - gamma_neg = 0.0_pReal + dot_gamma_neg = 0.0_pReal end where significantNegativeTau - if (present(dgamma_dtau_neg)) then + if (present(ddot_gamma_dtau_neg)) then significantNegativeTau2: where(abs(tau_neg)-dst%threshold_stress(:,of) > tol_math_check) dtn = t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_0 - dtk = t_k / tau_pos + dtk = t_k / tau_neg dvel = prm%kink_height * (dtk + dtn) / (t_n + t_k)**2.0_pReal - dgamma_dtau_neg = DotGamma0 * dvel * 0.5_pReal + ddot_gamma_dtau_neg = dot_gamma_0 * dvel * 0.5_pReal else where significantNegativeTau2 - dgamma_dtau_neg = 0.0_pReal + ddot_gamma_dtau_neg = 0.0_pReal end where significantNegativeTau2 end if From c9c9079076cce9385649fcf5e298889b0c3ae8cf Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 22 Mar 2019 10:55:17 +0100 Subject: [PATCH 10/12] fixed indentation --- src/plastic_disloUCLA.f90 | 1186 ++++++++++++++++++------------------- 1 file changed, 593 insertions(+), 593 deletions(-) diff --git a/src/plastic_disloUCLA.f90 b/src/plastic_disloUCLA.f90 index 22bb8ff96..e3b4bea18 100644 --- a/src/plastic_disloUCLA.f90 +++ b/src/plastic_disloUCLA.f90 @@ -6,102 +6,102 @@ !> @brief crystal plasticity model for bcc metals, especially Tungsten !-------------------------------------------------------------------------------------------------- module plastic_disloUCLA - use prec, only: & - pReal - - implicit none - private - integer, dimension(:,:), allocatable, target, public :: & - plastic_disloUCLA_sizePostResult !< size of each post result output - character(len=64), dimension(:,:), allocatable, target, public :: & - plastic_disloUCLA_output !< name of each post result output - - real(pReal), parameter, private :: & - kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin - - enum, bind(c) - enumerator :: & - undefined_ID, & - rho_mob_ID, & - rho_dip_ID, & - dot_gamma_sl_ID, & - gamma_sl_ID, & - Lambda_sl_ID, & - thresholdstress_ID - end enum - - type, private :: tParameters - real(pReal) :: & - aTol_rho, & - D, & !< grain size - mu, & - D_0, & !< prefactor for self-diffusion coefficient - Q_cl !< activation energy for dislocation climb - real(pReal), dimension(:), allocatable :: & - rho_mob_0, & !< initial dislocation density - rho_dip_0, & !< initial dipole density - b_sl, & !< magnitude of burgers vector [m] - nonSchmidCoeff, & - D_a, & - i_sl, & !< Adj. parameter for distance between 2 forest dislocations - atomicVolume, & - tau_0, & - !* mobility law parameters - delta_F, & !< activation energy for glide [J] - v0, & !< dislocation velocity prefactor [m/s] - p, & !< p-exponent in glide velocity - q, & !< q-exponent in glide velocity - B, & !< friction coefficient - kink_height, & !< height of the kink pair - w, & !< width of the kink pair - omega !< attempt frequency for kink pair nucleation - real(pReal), dimension(:,:), allocatable :: & - h_sl_sl, & !< slip resistance from slip activity - forestProjectionEdge - real(pReal), dimension(:,:,:), allocatable :: & - Schmid, & - nonSchmid_pos, & - nonSchmid_neg - integer :: & - sum_N_sl !< total number of active slip system - integer, dimension(:), allocatable :: & - N_sl !< number of active slip systems for each family - integer(kind(undefined_ID)), dimension(:),allocatable :: & - outputID !< ID of each post result output - logical :: & - dipoleFormation !< flag indicating consideration of dipole formation - end type !< container type for internal constitutive parameters - - type, private :: tDisloUCLAState - real(pReal), dimension(:,:), pointer :: & - rho_mob, & - rho_dip, & - gamma_sl - end type tDisloUCLAState - - type, private :: tDisloUCLAdependentState - real(pReal), dimension(:,:), allocatable :: & - Lambda_sl, & - threshold_stress - end type tDisloUCLAdependentState + use prec, only: & + pReal + + implicit none + private + integer, dimension(:,:), allocatable, target, public :: & + plastic_disloUCLA_sizePostResult !< size of each post result output + character(len=64), dimension(:,:), allocatable, target, public :: & + plastic_disloUCLA_output !< name of each post result output + + real(pReal), parameter, private :: & + kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin + + enum, bind(c) + enumerator :: & + undefined_ID, & + rho_mob_ID, & + rho_dip_ID, & + dot_gamma_sl_ID, & + gamma_sl_ID, & + Lambda_sl_ID, & + thresholdstress_ID + end enum + + type, private :: tParameters + real(pReal) :: & + aTol_rho, & + D, & !< grain size + mu, & + D_0, & !< prefactor for self-diffusion coefficient + Q_cl !< activation energy for dislocation climb + real(pReal), dimension(:), allocatable :: & + rho_mob_0, & !< initial dislocation density + rho_dip_0, & !< initial dipole density + b_sl, & !< magnitude of burgers vector [m] + nonSchmidCoeff, & + D_a, & + i_sl, & !< Adj. parameter for distance between 2 forest dislocations + atomicVolume, & + tau_0, & + !* mobility law parameters + delta_F, & !< activation energy for glide [J] + v0, & !< dislocation velocity prefactor [m/s] + p, & !< p-exponent in glide velocity + q, & !< q-exponent in glide velocity + B, & !< friction coefficient + kink_height, & !< height of the kink pair + w, & !< width of the kink pair + omega !< attempt frequency for kink pair nucleation + real(pReal), dimension(:,:), allocatable :: & + h_sl_sl, & !< slip resistance from slip activity + forestProjectionEdge + real(pReal), dimension(:,:,:), allocatable :: & + Schmid, & + nonSchmid_pos, & + nonSchmid_neg + integer :: & + sum_N_sl !< total number of active slip system + integer, dimension(:), allocatable :: & + N_sl !< number of active slip systems for each family + integer(kind(undefined_ID)), dimension(:),allocatable :: & + outputID !< ID of each post result output + logical :: & + dipoleFormation !< flag indicating consideration of dipole formation + end type !< container type for internal constitutive parameters + + type, private :: tDisloUCLAState + real(pReal), dimension(:,:), pointer :: & + rho_mob, & + rho_dip, & + gamma_sl + end type tDisloUCLAState + + type, private :: tDisloUCLAdependentState + real(pReal), dimension(:,:), allocatable :: & + Lambda_sl, & + threshold_stress + end type tDisloUCLAdependentState !-------------------------------------------------------------------------------------------------- ! containers for parameters and state - type(tParameters), allocatable, dimension(:), private :: param - type(tDisloUCLAState), allocatable, dimension(:), private :: & - dotState, & - state - type(tDisloUCLAdependentState), allocatable, dimension(:), private :: dependentState - - public :: & - plastic_disloUCLA_init, & - plastic_disloUCLA_dependentState, & - plastic_disloUCLA_LpAndItsTangent, & - plastic_disloUCLA_dotState, & - plastic_disloUCLA_postResults, & - plastic_disloUCLA_results - private :: & - kinetics + type(tParameters), allocatable, dimension(:), private :: param + type(tDisloUCLAState), allocatable, dimension(:), private :: & + dotState, & + state + type(tDisloUCLAdependentState), allocatable, dimension(:), private :: dependentState + + public :: & + plastic_disloUCLA_init, & + plastic_disloUCLA_dependentState, & + plastic_disloUCLA_LpAndItsTangent, & + plastic_disloUCLA_dotState, & + plastic_disloUCLA_postResults, & + plastic_disloUCLA_results + private :: & + kinetics contains @@ -111,247 +111,247 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine plastic_disloUCLA_init() - use prec, only: & - pStringLen - use debug, only: & - debug_level,& - debug_constitutive,& - debug_levelBasic - use math, only: & - math_expand - use IO, only: & - IO_error - use material, only: & - phase_plasticity, & - phase_plasticityInstance, & - phase_Noutput, & - material_allocatePlasticState, & - PLASTICITY_DISLOUCLA_label, & - PLASTICITY_DISLOUCLA_ID, & - material_phase, & - plasticState - use config, only: & - config_phase - use lattice - - implicit none - integer :: & - Ninstance, & - p, i, & - NipcMyPhase, & - sizeState, sizeDotState, & - startIndex, endIndex - - integer, dimension(0), parameter :: emptyIntArray = [integer::] - real(pReal), dimension(0), parameter :: emptyRealArray = [real(pReal)::] - character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] - - integer(kind(undefined_ID)) :: & - outputID - - character(len=pStringLen) :: & - extmsg = '' - character(len=65536), dimension(:), allocatable :: & - outputs - - write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_DISLOUCLA_label//' init -+>>>' - - write(6,'(/,a)') ' Cereceda et al., International Journal of Plasticity 78:242–256, 2016' - write(6,'(a)') ' https://dx.doi.org/10.1016/j.ijplas.2015.09.002' - - Ninstance = count(phase_plasticity == PLASTICITY_DISLOUCLA_ID) - if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & - write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - - allocate(plastic_disloUCLA_sizePostResult(maxval(phase_Noutput),Ninstance),source=0) - allocate(plastic_disloUCLA_output(maxval(phase_Noutput),Ninstance)) - plastic_disloUCLA_output = '' - - allocate(param(Ninstance)) - allocate(state(Ninstance)) - allocate(dotState(Ninstance)) - allocate(dependentState(Ninstance)) - - - do p = 1, size(phase_plasticity) - if (phase_plasticity(p) /= PLASTICITY_DISLOUCLA_ID) cycle - associate(prm => param(phase_plasticityInstance(p)), & - dot => dotState(phase_plasticityInstance(p)), & - stt => state(phase_plasticityInstance(p)), & - dst => dependentState(phase_plasticityInstance(p)), & - config => config_phase(p)) + use prec, only: & + pStringLen + use debug, only: & + debug_level,& + debug_constitutive,& + debug_levelBasic + use math, only: & + math_expand + use IO, only: & + IO_error + use material, only: & + phase_plasticity, & + phase_plasticityInstance, & + phase_Noutput, & + material_allocatePlasticState, & + PLASTICITY_DISLOUCLA_label, & + PLASTICITY_DISLOUCLA_ID, & + material_phase, & + plasticState + use config, only: & + config_phase + use lattice + + implicit none + integer :: & + Ninstance, & + p, i, & + NipcMyPhase, & + sizeState, sizeDotState, & + startIndex, endIndex + + integer, dimension(0), parameter :: emptyIntArray = [integer::] + real(pReal), dimension(0), parameter :: emptyRealArray = [real(pReal)::] + character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] + + integer(kind(undefined_ID)) :: & + outputID + + character(len=pStringLen) :: & + extmsg = '' + character(len=65536), dimension(:), allocatable :: & + outputs + + write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_DISLOUCLA_label//' init -+>>>' + + write(6,'(/,a)') ' Cereceda et al., International Journal of Plasticity 78:242–256, 2016' + write(6,'(a)') ' https://dx.doi.org/10.1016/j.ijplas.2015.09.002' + + Ninstance = count(phase_plasticity == PLASTICITY_DISLOUCLA_ID) + if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & + write(6,'(a16,1x,i5,/)') '# instances:',Ninstance + + allocate(plastic_disloUCLA_sizePostResult(maxval(phase_Noutput),Ninstance),source=0) + allocate(plastic_disloUCLA_output(maxval(phase_Noutput),Ninstance)) + plastic_disloUCLA_output = '' + + allocate(param(Ninstance)) + allocate(state(Ninstance)) + allocate(dotState(Ninstance)) + allocate(dependentState(Ninstance)) + + + do p = 1, size(phase_plasticity) + if (phase_plasticity(p) /= PLASTICITY_DISLOUCLA_ID) cycle + associate(prm => param(phase_plasticityInstance(p)), & + dot => dotState(phase_plasticityInstance(p)), & + stt => state(phase_plasticityInstance(p)), & + dst => dependentState(phase_plasticityInstance(p)), & + config => config_phase(p)) !-------------------------------------------------------------------------------------------------- ! optional parameters that need to be defined - prm%mu = lattice_mu(p) - - prm%aTol_rho = config%getFloat('atol_rho') - - ! sanity checks - if (prm%aTol_rho <= 0.0_pReal) extmsg = trim(extmsg)//' atol_rho' + prm%mu = lattice_mu(p) + + prm%aTol_rho = config%getFloat('atol_rho') + + ! sanity checks + if (prm%aTol_rho <= 0.0_pReal) extmsg = trim(extmsg)//' atol_rho' !-------------------------------------------------------------------------------------------------- ! slip related parameters - prm%N_sl = config%getInts('nslip',defaultVal=emptyIntArray) - prm%sum_N_sl = sum(prm%N_sl) - slipActive: if (prm%sum_N_sl > 0) then - prm%Schmid = lattice_SchmidMatrix_slip(prm%N_sl,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%N_sl,prm%nonSchmidCoeff,+1) - prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%N_sl,prm%nonSchmidCoeff,-1) - else - prm%nonSchmid_pos = prm%Schmid - prm%nonSchmid_neg = prm%Schmid - endif - - prm%h_sl_sl = lattice_interaction_SlipBySlip(prm%N_sl, & - config%getFloats('interaction_slipslip'), & - config%getString('lattice_structure')) - prm%forestProjectionEdge = lattice_forestProjection(prm%N_sl,config%getString('lattice_structure'),& - config%getFloat('c/a',defaultVal=0.0_pReal)) - - prm%rho_mob_0 = config%getFloats('rhoedge0', requiredSize=size(prm%N_sl)) - prm%rho_dip_0 = config%getFloats('rhoedgedip0', requiredSize=size(prm%N_sl)) - prm%v0 = config%getFloats('v0', requiredSize=size(prm%N_sl)) - prm%b_sl = config%getFloats('slipburgers', requiredSize=size(prm%N_sl)) - prm%delta_F = config%getFloats('qedge', requiredSize=size(prm%N_sl)) - - prm%i_sl = config%getFloats('clambdaslip', requiredSize=size(prm%N_sl)) - prm%tau_0 = config%getFloats('tau_peierls', requiredSize=size(prm%N_sl)) - prm%p = config%getFloats('p_slip', requiredSize=size(prm%N_sl), & - defaultVal=[(1.0_pReal,i=1,size(prm%N_sl))]) - prm%q = config%getFloats('q_slip', requiredSize=size(prm%N_sl), & - defaultVal=[(1.0_pReal,i=1,size(prm%N_sl))]) - prm%kink_height = config%getFloats('kink_height', requiredSize=size(prm%N_sl)) - prm%w = config%getFloats('kink_width', requiredSize=size(prm%N_sl)) - prm%omega = config%getFloats('omega', requiredSize=size(prm%N_sl)) - prm%B = config%getFloats('friction_coeff', requiredSize=size(prm%N_sl)) - - prm%D = config%getFloat('grainsize') - prm%D_0 = config%getFloat('d0') - prm%Q_cl = config%getFloat('qsd') - prm%atomicVolume = config%getFloat('catomicvolume') * prm%b_sl**3.0_pReal - prm%D_a = config%getFloat('cedgedipmindistance') * prm%b_sl - prm%dipoleformation = config%getFloat('dipoleformationfactor') > 0.0_pReal !should be on by default, ToDo: change to /key/-type key - - ! expand: family => system - prm%rho_mob_0 = math_expand(prm%rho_mob_0, prm%N_sl) - prm%rho_dip_0 = math_expand(prm%rho_dip_0, prm%N_sl) - prm%q = math_expand(prm%q, prm%N_sl) - prm%p = math_expand(prm%p, prm%N_sl) - prm%delta_F = math_expand(prm%delta_F, prm%N_sl) - prm%b_sl = math_expand(prm%b_sl, prm%N_sl) - prm%kink_height = math_expand(prm%kink_height, prm%N_sl) - prm%w = math_expand(prm%w, prm%N_sl) - prm%omega = math_expand(prm%omega, prm%N_sl) - prm%tau_0 = math_expand(prm%tau_0, prm%N_sl) - prm%v0 = math_expand(prm%v0, prm%N_sl) - prm%B = math_expand(prm%B, prm%N_sl) - prm%i_sl = math_expand(prm%i_sl, prm%N_sl) - prm%atomicVolume = math_expand(prm%atomicVolume, prm%N_sl) - prm%D_a = math_expand(prm%D_a, prm%N_sl) - - - ! sanity checks - if ( prm%D_0 <= 0.0_pReal) extmsg = trim(extmsg)//' D_0' - if ( prm%Q_cl <= 0.0_pReal) extmsg = trim(extmsg)//' Q_cl' - if (any(prm%rho_mob_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoedge0' - if (any(prm%rho_dip_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoedgedip0' - if (any(prm%v0 < 0.0_pReal)) extmsg = trim(extmsg)//' v0' - if (any(prm%b_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' slipb_sl' - if (any(prm%delta_F <= 0.0_pReal)) extmsg = trim(extmsg)//' qedge' - if (any(prm%tau_0 < 0.0_pReal)) extmsg = trim(extmsg)//' tau_0' - if (any(prm%D_a <= 0.0_pReal)) extmsg = trim(extmsg)//' cedgedipmindistance or slipb_sl' - if (any(prm%atomicVolume <= 0.0_pReal)) extmsg = trim(extmsg)//' catomicvolume or slipb_sl' - - else slipActive - allocate(prm%rho_mob_0(0)) - allocate(prm%rho_dip_0(0)) - endif slipActive + prm%N_sl = config%getInts('nslip',defaultVal=emptyIntArray) + prm%sum_N_sl = sum(prm%N_sl) + slipActive: if (prm%sum_N_sl > 0) then + prm%Schmid = lattice_SchmidMatrix_slip(prm%N_sl,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%N_sl,prm%nonSchmidCoeff,+1) + prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%N_sl,prm%nonSchmidCoeff,-1) + else + prm%nonSchmid_pos = prm%Schmid + prm%nonSchmid_neg = prm%Schmid + endif + + prm%h_sl_sl = lattice_interaction_SlipBySlip(prm%N_sl, & + config%getFloats('interaction_slipslip'), & + config%getString('lattice_structure')) + prm%forestProjectionEdge = lattice_forestProjection(prm%N_sl,config%getString('lattice_structure'),& + config%getFloat('c/a',defaultVal=0.0_pReal)) + + prm%rho_mob_0 = config%getFloats('rhoedge0', requiredSize=size(prm%N_sl)) + prm%rho_dip_0 = config%getFloats('rhoedgedip0', requiredSize=size(prm%N_sl)) + prm%v0 = config%getFloats('v0', requiredSize=size(prm%N_sl)) + prm%b_sl = config%getFloats('slipburgers', requiredSize=size(prm%N_sl)) + prm%delta_F = config%getFloats('qedge', requiredSize=size(prm%N_sl)) + + prm%i_sl = config%getFloats('clambdaslip', requiredSize=size(prm%N_sl)) + prm%tau_0 = config%getFloats('tau_peierls', requiredSize=size(prm%N_sl)) + prm%p = config%getFloats('p_slip', requiredSize=size(prm%N_sl), & + defaultVal=[(1.0_pReal,i=1,size(prm%N_sl))]) + prm%q = config%getFloats('q_slip', requiredSize=size(prm%N_sl), & + defaultVal=[(1.0_pReal,i=1,size(prm%N_sl))]) + prm%kink_height = config%getFloats('kink_height', requiredSize=size(prm%N_sl)) + prm%w = config%getFloats('kink_width', requiredSize=size(prm%N_sl)) + prm%omega = config%getFloats('omega', requiredSize=size(prm%N_sl)) + prm%B = config%getFloats('friction_coeff', requiredSize=size(prm%N_sl)) + + prm%D = config%getFloat('grainsize') + prm%D_0 = config%getFloat('d0') + prm%Q_cl = config%getFloat('qsd') + prm%atomicVolume = config%getFloat('catomicvolume') * prm%b_sl**3.0_pReal + prm%D_a = config%getFloat('cedgedipmindistance') * prm%b_sl + prm%dipoleformation = config%getFloat('dipoleformationfactor') > 0.0_pReal !should be on by default, ToDo: change to /key/-type key + + ! expand: family => system + prm%rho_mob_0 = math_expand(prm%rho_mob_0, prm%N_sl) + prm%rho_dip_0 = math_expand(prm%rho_dip_0, prm%N_sl) + prm%q = math_expand(prm%q, prm%N_sl) + prm%p = math_expand(prm%p, prm%N_sl) + prm%delta_F = math_expand(prm%delta_F, prm%N_sl) + prm%b_sl = math_expand(prm%b_sl, prm%N_sl) + prm%kink_height = math_expand(prm%kink_height, prm%N_sl) + prm%w = math_expand(prm%w, prm%N_sl) + prm%omega = math_expand(prm%omega, prm%N_sl) + prm%tau_0 = math_expand(prm%tau_0, prm%N_sl) + prm%v0 = math_expand(prm%v0, prm%N_sl) + prm%B = math_expand(prm%B, prm%N_sl) + prm%i_sl = math_expand(prm%i_sl, prm%N_sl) + prm%atomicVolume = math_expand(prm%atomicVolume, prm%N_sl) + prm%D_a = math_expand(prm%D_a, prm%N_sl) + + + ! sanity checks + if ( prm%D_0 <= 0.0_pReal) extmsg = trim(extmsg)//' D_0' + if ( prm%Q_cl <= 0.0_pReal) extmsg = trim(extmsg)//' Q_cl' + if (any(prm%rho_mob_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoedge0' + if (any(prm%rho_dip_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoedgedip0' + if (any(prm%v0 < 0.0_pReal)) extmsg = trim(extmsg)//' v0' + if (any(prm%b_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' slipb_sl' + if (any(prm%delta_F <= 0.0_pReal)) extmsg = trim(extmsg)//' qedge' + if (any(prm%tau_0 < 0.0_pReal)) extmsg = trim(extmsg)//' tau_0' + if (any(prm%D_a <= 0.0_pReal)) extmsg = trim(extmsg)//' cedgedipmindistance or slipb_sl' + if (any(prm%atomicVolume <= 0.0_pReal)) extmsg = trim(extmsg)//' catomicvolume or slipb_sl' + + else slipActive + allocate(prm%rho_mob_0(0)) + allocate(prm%rho_dip_0(0)) + endif slipActive !-------------------------------------------------------------------------------------------------- ! exit if any parameter is out of range - if (extmsg /= '') & - call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_DISLOUCLA_label//')') + if (extmsg /= '') & + call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_DISLOUCLA_label//')') !-------------------------------------------------------------------------------------------------- ! output pararameters - outputs = config%getStrings('(output)',defaultVal=emptyStringArray) - allocate(prm%outputID(0)) - do i=1, size(outputs) - outputID = undefined_ID - select case(trim(outputs(i))) - - case ('edge_density') - outputID = merge(rho_mob_ID,undefined_ID,prm%sum_N_sl>0) - case ('dipole_density') - outputID = merge(rho_dip_ID,undefined_ID,prm%sum_N_sl>0) - case ('shear_rate','shearrate','shear_rate_slip','shearrate_slip') - outputID = merge(dot_gamma_sl_ID,undefined_ID,prm%sum_N_sl>0) - case ('accumulated_shear','accumulatedshear','accumulated_shear_slip') - outputID = merge(gamma_sl_ID,undefined_ID,prm%sum_N_sl>0) - case ('mfp','mfp_slip') - outputID = merge(Lambda_sl_ID,undefined_ID,prm%sum_N_sl>0) - case ('threshold_stress','threshold_stress_slip') - outputID = merge(thresholdstress_ID,undefined_ID,prm%sum_N_sl>0) - - end select - - if (outputID /= undefined_ID) then - plastic_disloUCLA_output(i,phase_plasticityInstance(p)) = outputs(i) - plastic_disloUCLA_sizePostResult(i,phase_plasticityInstance(p)) = prm%sum_N_sl - prm%outputID = [prm%outputID, outputID] - endif - - enddo + outputs = config%getStrings('(output)',defaultVal=emptyStringArray) + allocate(prm%outputID(0)) + do i=1, size(outputs) + outputID = undefined_ID + select case(trim(outputs(i))) + + case ('edge_density') + outputID = merge(rho_mob_ID,undefined_ID,prm%sum_N_sl>0) + case ('dipole_density') + outputID = merge(rho_dip_ID,undefined_ID,prm%sum_N_sl>0) + case ('shear_rate','shearrate','shear_rate_slip','shearrate_slip') + outputID = merge(dot_gamma_sl_ID,undefined_ID,prm%sum_N_sl>0) + case ('accumulated_shear','accumulatedshear','accumulated_shear_slip') + outputID = merge(gamma_sl_ID,undefined_ID,prm%sum_N_sl>0) + case ('mfp','mfp_slip') + outputID = merge(Lambda_sl_ID,undefined_ID,prm%sum_N_sl>0) + case ('threshold_stress','threshold_stress_slip') + outputID = merge(thresholdstress_ID,undefined_ID,prm%sum_N_sl>0) + + end select + + if (outputID /= undefined_ID) then + plastic_disloUCLA_output(i,phase_plasticityInstance(p)) = outputs(i) + plastic_disloUCLA_sizePostResult(i,phase_plasticityInstance(p)) = prm%sum_N_sl + prm%outputID = [prm%outputID, outputID] + endif + + enddo !-------------------------------------------------------------------------------------------------- ! allocate state arrays - NipcMyPhase = count(material_phase == p) - sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl - sizeState = sizeDotState - - call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0, & - prm%sum_N_sl,0,0) - plasticState(p)%sizePostResults = sum(plastic_disloUCLA_sizePostResult(:,phase_plasticityInstance(p))) + NipcMyPhase = count(material_phase == p) + sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl + sizeState = sizeDotState + + call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0, & + prm%sum_N_sl,0,0) + plasticState(p)%sizePostResults = sum(plastic_disloUCLA_sizePostResult(:,phase_plasticityInstance(p))) !-------------------------------------------------------------------------------------------------- ! locally defined state aliases and initialization of state0 and aTolState - startIndex = 1 - endIndex = prm%sum_N_sl - stt%rho_mob=>plasticState(p)%state(startIndex:endIndex,:) - stt%rho_mob= spread(prm%rho_mob_0,2,NipcMyPhase) - dot%rho_mob=>plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%aTolState(startIndex:endIndex) = prm%aTol_rho - - startIndex = endIndex + 1 - endIndex = endIndex + prm%sum_N_sl - stt%rho_dip=>plasticState(p)%state(startIndex:endIndex,:) - stt%rho_dip= spread(prm%rho_dip_0,2,NipcMyPhase) - dot%rho_dip=>plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%aTolState(startIndex:endIndex) = prm%aTol_rho - - startIndex = endIndex + 1 - endIndex = endIndex + prm%sum_N_sl - stt%gamma_sl=>plasticState(p)%state(startIndex:endIndex,:) - dot%gamma_sl=>plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%aTolState(startIndex:endIndex) = 1.0e6_pReal !Don't use for convergence check - ! global alias - plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%accumulatedSlip => plasticState(p)%state(startIndex:endIndex,:) - - allocate(dst%Lambda_sl(prm%sum_N_sl,NipcMyPhase), source=0.0_pReal) - allocate(dst%threshold_stress(prm%sum_N_sl,NipcMyPhase), source=0.0_pReal) - - plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally - - end associate - - enddo + startIndex = 1 + endIndex = prm%sum_N_sl + stt%rho_mob=>plasticState(p)%state(startIndex:endIndex,:) + stt%rho_mob= spread(prm%rho_mob_0,2,NipcMyPhase) + dot%rho_mob=>plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%aTolState(startIndex:endIndex) = prm%aTol_rho + + startIndex = endIndex + 1 + endIndex = endIndex + prm%sum_N_sl + stt%rho_dip=>plasticState(p)%state(startIndex:endIndex,:) + stt%rho_dip= spread(prm%rho_dip_0,2,NipcMyPhase) + dot%rho_dip=>plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%aTolState(startIndex:endIndex) = prm%aTol_rho + + startIndex = endIndex + 1 + endIndex = endIndex + prm%sum_N_sl + stt%gamma_sl=>plasticState(p)%state(startIndex:endIndex,:) + dot%gamma_sl=>plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%aTolState(startIndex:endIndex) = 1.0e6_pReal ! Don't use for convergence check + ! global alias + plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%accumulatedSlip => plasticState(p)%state(startIndex:endIndex,:) + + allocate(dst%Lambda_sl(prm%sum_N_sl,NipcMyPhase), source=0.0_pReal) + allocate(dst%threshold_stress(prm%sum_N_sl,NipcMyPhase), source=0.0_pReal) + + plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally + + end associate + + enddo end subroutine plastic_disloUCLA_init @@ -359,43 +359,43 @@ end subroutine plastic_disloUCLA_init !-------------------------------------------------------------------------------------------------- !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- -pure subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) - - implicit none - real(pReal), dimension(3,3), intent(out) :: & - Lp !< plastic velocity gradient - real(pReal), dimension(3,3,3,3), intent(out) :: & - dLp_dMp !< derivative of Lp with respect to the Mandel stress - - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - real(pReal), intent(in) :: & - T !< temperature - integer, intent(in) :: & - instance, & - of - - integer :: & - i,k,l,m,n - real(pReal), dimension(param(instance)%sum_N_sl) :: & - dot_gamma_pos,dot_gamma_neg, & - ddot_gamma_dtau_pos,ddot_gamma_dtau_neg - - Lp = 0.0_pReal - dLp_dMp = 0.0_pReal - - associate(prm => param(instance)) - - call kinetics(Mp,T,instance,of,dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg) - do i = 1, prm%sum_N_sl - Lp = Lp + (dot_gamma_pos(i)+dot_gamma_neg(i))*prm%Schmid(1:3,1:3,i) - forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & - + ddot_gamma_dtau_pos(i) * prm%Schmid(k,l,i) * prm%nonSchmid_pos(m,n,i) & - + ddot_gamma_dtau_neg(i) * prm%Schmid(k,l,i) * prm%nonSchmid_neg(m,n,i) - enddo - - end associate +pure subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp, & + Mp,T,instance,of) + implicit none + real(pReal), dimension(3,3), intent(out) :: & + Lp !< plastic velocity gradient + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp !< derivative of Lp with respect to the Mandel stress + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + T !< temperature + integer, intent(in) :: & + instance, & + of + + integer :: & + i,k,l,m,n + real(pReal), dimension(param(instance)%sum_N_sl) :: & + dot_gamma_pos,dot_gamma_neg, & + ddot_gamma_dtau_pos,ddot_gamma_dtau_neg + + Lp = 0.0_pReal + dLp_dMp = 0.0_pReal + + associate(prm => param(instance)) + + call kinetics(Mp,T,instance,of,dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg) + do i = 1, prm%sum_N_sl + Lp = Lp + (dot_gamma_pos(i)+dot_gamma_neg(i))*prm%Schmid(1:3,1:3,i) + forall (k=1:3,l=1:3,m=1:3,n=1:3) & + dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + + ddot_gamma_dtau_pos(i) * prm%Schmid(k,l,i) * prm%nonSchmid_pos(m,n,i) & + + ddot_gamma_dtau_neg(i) * prm%Schmid(k,l,i) * prm%nonSchmid_neg(m,n,i) + enddo + + end associate end subroutine plastic_disloUCLA_LpAndItsTangent @@ -404,63 +404,63 @@ end subroutine plastic_disloUCLA_LpAndItsTangent !> @brief calculates the rate of change of microstructure !-------------------------------------------------------------------------------------------------- subroutine plastic_disloUCLA_dotState(Mp,T,instance,of) - use prec, only: & - tol_math_check, & - dEq0 - use math, only: & - PI, & - math_clip - - implicit none - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - real(pReal), intent(in) :: & - T !< temperature - integer, intent(in) :: & - instance, & - of - - real(pReal) :: & - VacancyDiffusion - real(pReal), dimension(param(instance)%sum_N_sl) :: & - gdot_pos, gdot_neg,& - tau_pos,& - tau_neg, & - DotRhoDipFormation, v_cl, EdgeDipDistance, & - DotRhoEdgeDipClimb - - associate(prm => param(instance), stt => state(instance),dot => dotState(instance), dst => dependentState(instance)) - - call kinetics(Mp,T,instance,of,& - gdot_pos,gdot_neg, & - tau_pos_out = tau_pos,tau_neg_out = tau_neg) - - dot%gamma_sl(:,of) = (gdot_pos+gdot_neg) ! ToDo: needs to be abs - VacancyDiffusion = prm%D_0*exp(-prm%Q_cl/(kB*T)) - - where(dEq0(tau_pos)) ! ToDo: use avg of pos and neg - DotRhoDipFormation = 0.0_pReal - DotRhoEdgeDipClimb = 0.0_pReal - else where - EdgeDipDistance = math_clip((3.0_pReal*prm%mu*prm%b_sl)/(16.0_pReal*PI*abs(tau_pos)), & - prm%D_a, & ! lower limit - dst%Lambda_sl(:,of)) ! upper limit - dotRhoDipFormation = merge(((2.0_pReal*EdgeDipDistance)/prm%b_sl)* stt%rho_mob(:,of)*abs(dot%gamma_sl(:,of)), & ! ToDo: ignore region of spontaneous annihilation - 0.0_pReal, & - prm%dipoleformation) - v_cl = (3.0_pReal*prm%mu*VacancyDiffusion*prm%atomicVolume/(2.0_pReal*pi*kB*T)) & - * (1.0_pReal/(EdgeDipDistance+prm%D_a)) - dotRhoEdgeDipClimb = (4.0_pReal*v_cl*stt%rho_dip(:,of))/(EdgeDipDistance-prm%D_a) ! ToDo: Discuss with Franz: Stress dependency? - end where - - dot%rho_mob(:,of) = abs(dot%gamma_sl(:,of))/(prm%b_sl*dst%Lambda_sl(:,of)) & ! multiplication - - DotRhoDipFormation & - - (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_mob(:,of)*abs(dot%gamma_sl(:,of)) ! Spontaneous annihilation of 2 single edge dislocations - dot%rho_dip(:,of) = DotRhoDipFormation & - - (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_dip(:,of)*abs(dot%gamma_sl(:,of)) & ! Spontaneous annihilation of a single edge dislocation with a dipole constituent - - DotRhoEdgeDipClimb - - end associate + use prec, only: & + tol_math_check, & + dEq0 + use math, only: & + PI, & + math_clip + + implicit none + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + T !< temperature + integer, intent(in) :: & + instance, & + of + + real(pReal) :: & + VacancyDiffusion + real(pReal), dimension(param(instance)%sum_N_sl) :: & + gdot_pos, gdot_neg,& + tau_pos,& + tau_neg, & + DotRhoDipFormation, v_cl, EdgeDipDistance, & + DotRhoEdgeDipClimb + + associate(prm => param(instance), stt => state(instance),dot => dotState(instance), dst => dependentState(instance)) + + call kinetics(Mp,T,instance,of,& + gdot_pos,gdot_neg, & + tau_pos_out = tau_pos,tau_neg_out = tau_neg) + + dot%gamma_sl(:,of) = (gdot_pos+gdot_neg) ! ToDo: needs to be abs + VacancyDiffusion = prm%D_0*exp(-prm%Q_cl/(kB*T)) + + where(dEq0(tau_pos)) ! ToDo: use avg of pos and neg + DotRhoDipFormation = 0.0_pReal + DotRhoEdgeDipClimb = 0.0_pReal + else where + EdgeDipDistance = math_clip((3.0_pReal*prm%mu*prm%b_sl)/(16.0_pReal*PI*abs(tau_pos)), & + prm%D_a, & ! lower limit + dst%Lambda_sl(:,of)) ! upper limit + dotRhoDipFormation = merge(((2.0_pReal*EdgeDipDistance)/prm%b_sl)* stt%rho_mob(:,of)*abs(dot%gamma_sl(:,of)), & ! ToDo: ignore region of spontaneous annihilation + 0.0_pReal, & + prm%dipoleformation) + v_cl = (3.0_pReal*prm%mu*VacancyDiffusion*prm%atomicVolume/(2.0_pReal*pi*kB*T)) & + * (1.0_pReal/(EdgeDipDistance+prm%D_a)) + dotRhoEdgeDipClimb = (4.0_pReal*v_cl*stt%rho_dip(:,of))/(EdgeDipDistance-prm%D_a) ! ToDo: Discuss with Franz: Stress dependency? + end where + + dot%rho_mob(:,of) = abs(dot%gamma_sl(:,of))/(prm%b_sl*dst%Lambda_sl(:,of)) & ! multiplication + - DotRhoDipFormation & + - (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_mob(:,of)*abs(dot%gamma_sl(:,of)) ! Spontaneous annihilation of 2 single edge dislocations + dot%rho_dip(:,of) = DotRhoDipFormation & + - (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_dip(:,of)*abs(dot%gamma_sl(:,of)) & ! Spontaneous annihilation of a single edge dislocation with a dipole constituent + - DotRhoEdgeDipClimb + + end associate end subroutine plastic_disloUCLA_dotState @@ -469,30 +469,30 @@ end subroutine plastic_disloUCLA_dotState !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- subroutine plastic_disloUCLA_dependentState(instance,of) - - implicit none - integer, intent(in) :: & - instance, & - of - - real(pReal), dimension(param(instance)%sum_N_sl) :: & - dislocationSpacing - integer :: & - i - - associate(prm => param(instance), stt => state(instance),dst => dependentState(instance)) - - forall (i = 1:prm%sum_N_sl) - dislocationSpacing(i) = sqrt(dot_product(stt%rho_mob(:,of)+stt%rho_dip(:,of), & - prm%forestProjectionEdge(:,i))) - dst%threshold_stress(i,of) = prm%mu*prm%b_sl(i) & - * sqrt(dot_product(stt%rho_mob(:,of)+stt%rho_dip(:,of), & - prm%h_sl_sl(:,i))) - end forall - - dst%Lambda_sl(:,of) = prm%D/(1.0_pReal+prm%D*dislocationSpacing/prm%i_sl) - - end associate + + implicit none + integer, intent(in) :: & + instance, & + of + + real(pReal), dimension(param(instance)%sum_N_sl) :: & + dislocationSpacing + integer :: & + i + + associate(prm => param(instance), stt => state(instance),dst => dependentState(instance)) + + forall (i = 1:prm%sum_N_sl) + dislocationSpacing(i) = sqrt(dot_product(stt%rho_mob(:,of)+stt%rho_dip(:,of), & + prm%forestProjectionEdge(:,i))) + dst%threshold_stress(i,of) = prm%mu*prm%b_sl(i) & + * sqrt(dot_product(stt%rho_mob(:,of)+stt%rho_dip(:,of), & + prm%h_sl_sl(:,i))) + end forall + + dst%Lambda_sl(:,of) = prm%D/(1.0_pReal+prm%D*dislocationSpacing/prm%i_sl) + + end associate end subroutine plastic_disloUCLA_dependentState @@ -501,57 +501,57 @@ end subroutine plastic_disloUCLA_dependentState !> @brief return array of constitutive results !-------------------------------------------------------------------------------------------------- function plastic_disloUCLA_postResults(Mp,T,instance,of) result(postResults) - use prec, only: & - dEq, dNeq0 - use math, only: & - PI, & - math_mul33xx33 + use prec, only: & + dEq, dNeq0 + use math, only: & + PI, & + math_mul33xx33 + + implicit none + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + T !< temperature + integer, intent(in) :: & + instance, & + of + + real(pReal), dimension(sum(plastic_disloUCLA_sizePostResult(:,instance))) :: & + postResults + + integer :: & + o,c,i + real(pReal), dimension(param(instance)%sum_N_sl) :: & + gdot_pos,gdot_neg + + c = 0 + + associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) + + outputsLoop: do o = 1,size(prm%outputID) + select case(prm%outputID(o)) + + case (rho_mob_ID) + postResults(c+1:c+prm%sum_N_sl) = stt%rho_mob(1:prm%sum_N_sl,of) + case (rho_dip_ID) + postResults(c+1:c+prm%sum_N_sl) = stt%rho_dip(1:prm%sum_N_sl,of) + case (dot_gamma_sl_ID) + call kinetics(Mp,T,instance,of,gdot_pos,gdot_neg) + postResults(c+1:c+prm%sum_N_sl) = gdot_pos + gdot_neg + case (gamma_sl_ID) + postResults(c+1:c+prm%sum_N_sl) = stt%gamma_sl(1:prm%sum_N_sl, of) + case (Lambda_sl_ID) + postResults(c+1:c+prm%sum_N_sl) = dst%Lambda_sl(1:prm%sum_N_sl, of) + case (thresholdstress_ID) + postResults(c+1:c+prm%sum_N_sl) = dst%threshold_stress(1:prm%sum_N_sl,of) + + end select + + c = c + prm%sum_N_sl + + enddo outputsLoop - implicit none - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - real(pReal), intent(in) :: & - T !< temperature - integer, intent(in) :: & - instance, & - of - - real(pReal), dimension(sum(plastic_disloUCLA_sizePostResult(:,instance))) :: & - postResults - - integer :: & - o,c,i - real(pReal), dimension(param(instance)%sum_N_sl) :: & - gdot_pos,gdot_neg - - c = 0 - - associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) - - outputsLoop: do o = 1,size(prm%outputID) - select case(prm%outputID(o)) - - case (rho_mob_ID) - postResults(c+1:c+prm%sum_N_sl) = stt%rho_mob(1:prm%sum_N_sl,of) - case (rho_dip_ID) - postResults(c+1:c+prm%sum_N_sl) = stt%rho_dip(1:prm%sum_N_sl,of) - case (dot_gamma_sl_ID) - call kinetics(Mp,T,instance,of,gdot_pos,gdot_neg) - postResults(c+1:c+prm%sum_N_sl) = gdot_pos + gdot_neg - case (gamma_sl_ID) - postResults(c+1:c+prm%sum_N_sl) = stt%gamma_sl(1:prm%sum_N_sl, of) - case (Lambda_sl_ID) - postResults(c+1:c+prm%sum_N_sl) = dst%Lambda_sl(1:prm%sum_N_sl, of) - case (thresholdstress_ID) - postResults(c+1:c+prm%sum_N_sl) = dst%threshold_stress(1:prm%sum_N_sl,of) - - end select - - c = c + prm%sum_N_sl - - enddo outputsLoop - - end associate + end associate end function plastic_disloUCLA_postResults @@ -591,116 +591,116 @@ end subroutine plastic_disloUCLA_results !-------------------------------------------------------------------------------------------------- pure subroutine kinetics(Mp,T,instance,of, & dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg,tau_pos_out,tau_neg_out) - use prec, only: & - tol_math_check, & - dEq, dNeq0 - use math, only: & - PI, & - math_mul33xx33 - - implicit none - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - real(pReal), intent(in) :: & - T !< temperature - integer, intent(in) :: & - instance, & - of - - real(pReal), intent(out), dimension(param(instance)%sum_N_sl) :: & - dot_gamma_pos, & - dot_gamma_neg - real(pReal), intent(out), optional, dimension(param(instance)%sum_N_sl) :: & - ddot_gamma_dtau_pos, & - ddot_gamma_dtau_neg, & - tau_pos_out, & - tau_neg_out - real(pReal), dimension(param(instance)%sum_N_sl) :: & - StressRatio, & - StressRatio_p,StressRatio_pminus1, & - dvel, vel, & - tau_pos,tau_neg, & - t_n, t_k, dtk,dtn, & - needsGoodName ! ToDo: @Karo: any idea? - integer :: j - - associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) - - do j = 1, prm%sum_N_sl - tau_pos(j) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,j)) - tau_neg(j) = math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,j)) - enddo - - - if (present(tau_pos_out)) tau_pos_out = tau_pos - if (present(tau_neg_out)) tau_neg_out = tau_neg - - associate(BoltzmannRatio => prm%delta_F/(kB*T), & - dot_gamma_0 => stt%rho_mob(:,of)*prm%b_sl*prm%v0, & - effectiveLength => dst%Lambda_sl(:,of) - prm%w) - - significantPositiveTau: where(abs(tau_pos)-dst%threshold_stress(:,of) > tol_math_check) - StressRatio = (abs(tau_pos)-dst%threshold_stress(:,of))/prm%tau_0 - StressRatio_p = StressRatio** prm%p - StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) - needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) - - t_n = prm%b_sl/(needsGoodName*prm%omega*effectiveLength) - t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos) ! our definition of tk is different with the one in dislotwin - - vel = prm%kink_height/(t_n + t_k) - - dot_gamma_pos = dot_gamma_0 * sign(vel,tau_pos) * 0.5_pReal - else where significantPositiveTau - dot_gamma_pos = 0.0_pReal - end where significantPositiveTau - - if (present(ddot_gamma_dtau_pos)) then - significantPositiveTau2: where(abs(tau_pos)-dst%threshold_stress(:,of) > tol_math_check) - dtn = t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & - * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_0 - dtk = t_k / tau_pos - - dvel = prm%kink_height * (dtk + dtn) / (t_n + t_k)**2.0_pReal - - ddot_gamma_dtau_pos = dot_gamma_0 * dvel* 0.5_pReal - else where significantPositiveTau2 - ddot_gamma_dtau_pos = 0.0_pReal - end where significantPositiveTau2 - endif - - significantNegativeTau: where(abs(tau_neg)-dst%threshold_stress(:,of) > tol_math_check) - StressRatio = (abs(tau_neg)-dst%threshold_stress(:,of))/prm%tau_0 - StressRatio_p = StressRatio** prm%p - StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) - needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) - - t_n = prm%b_sl/(needsGoodName*prm%omega*effectiveLength) - t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos) ! our definition of tk is different with the one in dislotwin - - vel = prm%kink_height/(t_n + t_k) - - dot_gamma_neg = dot_gamma_0 * sign(vel,tau_neg) * 0.5_pReal - else where significantNegativeTau - dot_gamma_neg = 0.0_pReal - end where significantNegativeTau - - if (present(ddot_gamma_dtau_neg)) then - significantNegativeTau2: where(abs(tau_neg)-dst%threshold_stress(:,of) > tol_math_check) - dtn = t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & - * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_0 - dtk = t_k / tau_neg - - dvel = prm%kink_height * (dtk + dtn) / (t_n + t_k)**2.0_pReal - - ddot_gamma_dtau_neg = dot_gamma_0 * dvel * 0.5_pReal - else where significantNegativeTau2 - ddot_gamma_dtau_neg = 0.0_pReal - end where significantNegativeTau2 - end if + use prec, only: & + tol_math_check, & + dEq, dNeq0 + use math, only: & + PI, & + math_mul33xx33 - end associate - end associate + implicit none + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + T !< temperature + integer, intent(in) :: & + instance, & + of + + real(pReal), intent(out), dimension(param(instance)%sum_N_sl) :: & + dot_gamma_pos, & + dot_gamma_neg + real(pReal), intent(out), optional, dimension(param(instance)%sum_N_sl) :: & + ddot_gamma_dtau_pos, & + ddot_gamma_dtau_neg, & + tau_pos_out, & + tau_neg_out + real(pReal), dimension(param(instance)%sum_N_sl) :: & + StressRatio, & + StressRatio_p,StressRatio_pminus1, & + dvel, vel, & + tau_pos,tau_neg, & + t_n, t_k, dtk,dtn, & + needsGoodName ! ToDo: @Karo: any idea? + integer :: j + + associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) + + do j = 1, prm%sum_N_sl + tau_pos(j) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,j)) + tau_neg(j) = math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,j)) + enddo + + + if (present(tau_pos_out)) tau_pos_out = tau_pos + if (present(tau_neg_out)) tau_neg_out = tau_neg + + associate(BoltzmannRatio => prm%delta_F/(kB*T), & + dot_gamma_0 => stt%rho_mob(:,of)*prm%b_sl*prm%v0, & + effectiveLength => dst%Lambda_sl(:,of) - prm%w) + + significantPositiveTau: where(abs(tau_pos)-dst%threshold_stress(:,of) > tol_math_check) + StressRatio = (abs(tau_pos)-dst%threshold_stress(:,of))/prm%tau_0 + StressRatio_p = StressRatio** prm%p + StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) + needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) + + t_n = prm%b_sl/(needsGoodName*prm%omega*effectiveLength) + t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos) ! our definition of tk is different with the one in dislotwin + + vel = prm%kink_height/(t_n + t_k) + + dot_gamma_pos = dot_gamma_0 * sign(vel,tau_pos) * 0.5_pReal + else where significantPositiveTau + dot_gamma_pos = 0.0_pReal + end where significantPositiveTau + + if (present(ddot_gamma_dtau_pos)) then + significantPositiveTau2: where(abs(tau_pos)-dst%threshold_stress(:,of) > tol_math_check) + dtn = t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & + * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_0 + dtk = t_k / tau_pos + + dvel = prm%kink_height * (dtk + dtn) / (t_n + t_k)**2.0_pReal + + ddot_gamma_dtau_pos = dot_gamma_0 * dvel* 0.5_pReal + else where significantPositiveTau2 + ddot_gamma_dtau_pos = 0.0_pReal + end where significantPositiveTau2 + endif + + significantNegativeTau: where(abs(tau_neg)-dst%threshold_stress(:,of) > tol_math_check) + StressRatio = (abs(tau_neg)-dst%threshold_stress(:,of))/prm%tau_0 + StressRatio_p = StressRatio** prm%p + StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) + needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) + + t_n = prm%b_sl/(needsGoodName*prm%omega*effectiveLength) + t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos) ! our definition of tk is different with the one in dislotwin + + vel = prm%kink_height/(t_n + t_k) + + dot_gamma_neg = dot_gamma_0 * sign(vel,tau_neg) * 0.5_pReal + else where significantNegativeTau + dot_gamma_neg = 0.0_pReal + end where significantNegativeTau + + if (present(ddot_gamma_dtau_neg)) then + significantNegativeTau2: where(abs(tau_neg)-dst%threshold_stress(:,of) > tol_math_check) + dtn = t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & + * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_0 + dtk = t_k / tau_neg + + dvel = prm%kink_height * (dtk + dtn) / (t_n + t_k)**2.0_pReal + + ddot_gamma_dtau_neg = dot_gamma_0 * dvel * 0.5_pReal + else where significantNegativeTau2 + ddot_gamma_dtau_neg = 0.0_pReal + end where significantNegativeTau2 + end if + + end associate + end associate end subroutine kinetics From bd6b7e185426f1ec9390048bca17fa19af83ed1b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 22 Mar 2019 10:58:54 +0100 Subject: [PATCH 11/12] correct signs (same behavior, but full math) --- src/plastic_disloUCLA.f90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/plastic_disloUCLA.f90 b/src/plastic_disloUCLA.f90 index e3b4bea18..6aebae603 100644 --- a/src/plastic_disloUCLA.f90 +++ b/src/plastic_disloUCLA.f90 @@ -657,11 +657,11 @@ pure subroutine kinetics(Mp,T,instance,of, & if (present(ddot_gamma_dtau_pos)) then significantPositiveTau2: where(abs(tau_pos)-dst%threshold_stress(:,of) > tol_math_check) - dtn = t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & + dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_0 - dtk = t_k / tau_pos + dtk = -1.0_pReal * t_k / tau_pos - dvel = prm%kink_height * (dtk + dtn) / (t_n + t_k)**2.0_pReal + dvel = -1.0_pReal * prm%kink_height * (dtk + dtn) / (t_n + t_k)**2.0_pReal ddot_gamma_dtau_pos = dot_gamma_0 * dvel* 0.5_pReal else where significantPositiveTau2 @@ -687,11 +687,11 @@ pure subroutine kinetics(Mp,T,instance,of, & if (present(ddot_gamma_dtau_neg)) then significantNegativeTau2: where(abs(tau_neg)-dst%threshold_stress(:,of) > tol_math_check) - dtn = t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & + dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_0 - dtk = t_k / tau_neg + dtk = -1.0_pReal * t_k / tau_neg - dvel = prm%kink_height * (dtk + dtn) / (t_n + t_k)**2.0_pReal + dvel = -1.0_pReal * prm%kink_height * (dtk + dtn) / (t_n + t_k)**2.0_pReal ddot_gamma_dtau_neg = dot_gamma_0 * dvel * 0.5_pReal else where significantNegativeTau2 From 450a0565d367f1b0fa1fa8650b209f311e3b9f03 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 22 Mar 2019 11:06:08 +0100 Subject: [PATCH 12/12] internal variable names also adjusted --- src/plastic_disloUCLA.f90 | 34 ++++++++++++++++++---------------- 1 file changed, 18 insertions(+), 16 deletions(-) diff --git a/src/plastic_disloUCLA.f90 b/src/plastic_disloUCLA.f90 index 6aebae603..1cdebe87f 100644 --- a/src/plastic_disloUCLA.f90 +++ b/src/plastic_disloUCLA.f90 @@ -426,8 +426,10 @@ subroutine plastic_disloUCLA_dotState(Mp,T,instance,of) gdot_pos, gdot_neg,& tau_pos,& tau_neg, & - DotRhoDipFormation, v_cl, EdgeDipDistance, & - DotRhoEdgeDipClimb + v_cl, & + dot_rho_dip_formation, & + dot_rho_dip_climb, & + dip_distance associate(prm => param(instance), stt => state(instance),dot => dotState(instance), dst => dependentState(instance)) @@ -439,26 +441,26 @@ subroutine plastic_disloUCLA_dotState(Mp,T,instance,of) VacancyDiffusion = prm%D_0*exp(-prm%Q_cl/(kB*T)) where(dEq0(tau_pos)) ! ToDo: use avg of pos and neg - DotRhoDipFormation = 0.0_pReal - DotRhoEdgeDipClimb = 0.0_pReal + dot_rho_dip_formation = 0.0_pReal + dot_rho_dip_climb = 0.0_pReal else where - EdgeDipDistance = math_clip((3.0_pReal*prm%mu*prm%b_sl)/(16.0_pReal*PI*abs(tau_pos)), & - prm%D_a, & ! lower limit - dst%Lambda_sl(:,of)) ! upper limit - dotRhoDipFormation = merge(((2.0_pReal*EdgeDipDistance)/prm%b_sl)* stt%rho_mob(:,of)*abs(dot%gamma_sl(:,of)), & ! ToDo: ignore region of spontaneous annihilation - 0.0_pReal, & - prm%dipoleformation) + dip_distance = math_clip(3.0_pReal*prm%mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos)), & + prm%D_a, & ! lower limit + dst%Lambda_sl(:,of)) ! upper limit + dot_rho_dip_formation = merge(2.0_pReal*dip_distance* stt%rho_mob(:,of)*abs(dot%gamma_sl(:,of))/prm%b_sl, & ! ToDo: ignore region of spontaneous annihilation + 0.0_pReal, & + prm%dipoleformation) v_cl = (3.0_pReal*prm%mu*VacancyDiffusion*prm%atomicVolume/(2.0_pReal*pi*kB*T)) & - * (1.0_pReal/(EdgeDipDistance+prm%D_a)) - dotRhoEdgeDipClimb = (4.0_pReal*v_cl*stt%rho_dip(:,of))/(EdgeDipDistance-prm%D_a) ! ToDo: Discuss with Franz: Stress dependency? + * (1.0_pReal/(dip_distance+prm%D_a)) + dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,of))/(dip_distance-prm%D_a) ! ToDo: Discuss with Franz: Stress dependency? end where dot%rho_mob(:,of) = abs(dot%gamma_sl(:,of))/(prm%b_sl*dst%Lambda_sl(:,of)) & ! multiplication - - DotRhoDipFormation & + - dot_rho_dip_formation & - (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_mob(:,of)*abs(dot%gamma_sl(:,of)) ! Spontaneous annihilation of 2 single edge dislocations - dot%rho_dip(:,of) = DotRhoDipFormation & + dot%rho_dip(:,of) = dot_rho_dip_formation & - (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_dip(:,of)*abs(dot%gamma_sl(:,of)) & ! Spontaneous annihilation of a single edge dislocation with a dipole constituent - - DotRhoEdgeDipClimb + - dot_rho_dip_climb end associate @@ -520,7 +522,7 @@ function plastic_disloUCLA_postResults(Mp,T,instance,of) result(postResults) postResults integer :: & - o,c,i + o,c real(pReal), dimension(param(instance)%sum_N_sl) :: & gdot_pos,gdot_neg