diff --git a/src/constitutive_plastic_disloUCLA.f90 b/src/constitutive_plastic_disloUCLA.f90 index 2b353ba26..00a1b8e21 100644 --- a/src/constitutive_plastic_disloUCLA.f90 +++ b/src/constitutive_plastic_disloUCLA.f90 @@ -6,21 +6,10 @@ !> @brief crystal plasticity model for bcc metals, especially Tungsten !-------------------------------------------------------------------------------------------------- submodule(constitutive) plastic_disloUCLA - + real(pReal), parameter :: & 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, & - tau_pass_ID - end enum - + type :: tParameters real(pReal) :: & aTol_rho, & @@ -28,7 +17,7 @@ submodule(constitutive) plastic_disloUCLA mu, & D_0, & !< prefactor for self-diffusion coefficient Q_cl !< activation energy for dislocation climb - real(pReal), dimension(:), allocatable :: & + real(pReal), allocatable, dimension(:) :: & rho_mob_0, & !< initial dislocation density rho_dip_0, & !< initial dipole density b_sl, & !< magnitude of burgers vector [m] @@ -46,30 +35,30 @@ submodule(constitutive) plastic_disloUCLA kink_height, & !< height of the kink pair w, & !< width of the kink pair omega !< attempt frequency for kink pair nucleation - real(pReal), dimension(:,:), allocatable :: & + real(pReal), allocatable, dimension(:,:) :: & h_sl_sl, & !< slip resistance from slip activity forestProjectionEdge - real(pReal), dimension(:,:,:), allocatable :: & + real(pReal), allocatable, dimension(:,:,:) :: & Schmid, & nonSchmid_pos, & nonSchmid_neg integer :: & sum_N_sl !< total number of active slip system - integer, dimension(:), allocatable :: & + integer, allocatable, dimension(:) :: & N_sl !< number of active slip systems for each family - integer(kind(undefined_ID)), dimension(:),allocatable :: & - outputID !< ID of each post result output + character(len=pStringLen), allocatable, dimension(:) :: & + output logical :: & dipoleFormation !< flag indicating consideration of dipole formation end type !< container type for internal constitutive parameters - + type :: tDisloUCLAState real(pReal), dimension(:,:), pointer :: & rho_mob, & rho_dip, & gamma_sl end type tDisloUCLAState - + type :: tDisloUCLAdependentState real(pReal), dimension(:,:), allocatable :: & Lambda_sl, & @@ -88,41 +77,36 @@ contains !-------------------------------------------------------------------------------------------------- -!> @brief module initialization +!> @brief Perform module initialization. !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- module subroutine plastic_disloUCLA_init - + integer :: & Ninstance, & p, i, & NipcMyPhase, & sizeState, sizeDotState, & startIndex, endIndex - - integer(kind(undefined_ID)) :: & - outputID - + character(len=pStringLen) :: & extmsg = '' - character(len=pStringLen), 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' - + + write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_DISLOUCLA_label//' init -+>>>'; flush(6) + + 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(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)), & @@ -134,9 +118,9 @@ module subroutine plastic_disloUCLA_init !-------------------------------------------------------------------------------------------------- ! 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' @@ -147,7 +131,7 @@ module subroutine plastic_disloUCLA_init 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) @@ -158,20 +142,20 @@ module subroutine plastic_disloUCLA_init 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_edge(prm%N_sl,config%getString('lattice_structure'),& config%getFloat('c/a',defaultVal=0.0_pReal)) prm%forestProjectionEdge = transpose(prm%forestProjectionEdge) - + 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), & @@ -182,14 +166,14 @@ module subroutine plastic_disloUCLA_init 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) @@ -206,8 +190,8 @@ module subroutine plastic_disloUCLA_init 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' @@ -219,7 +203,7 @@ module subroutine plastic_disloUCLA_init 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)) @@ -232,39 +216,14 @@ module subroutine plastic_disloUCLA_init !-------------------------------------------------------------------------------------------------- ! 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(tau_pass_ID,undefined_ID,prm%sum_N_sl>0) - - end select - - if (outputID /= undefined_ID) then - prm%outputID = [prm%outputID, outputID] - endif - - enddo + prm%output = config%getStrings('(output)',defaultVal=emptyStringArray) !-------------------------------------------------------------------------------------------------- ! allocate state arrays NipcMyPhase = count(material_phaseAt == p) * discretization_nIP sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl sizeState = sizeDotState - + call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0) !-------------------------------------------------------------------------------------------------- @@ -275,14 +234,14 @@ module subroutine plastic_disloUCLA_init 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,:) @@ -290,21 +249,21 @@ module subroutine plastic_disloUCLA_init plasticState(p)%aTolState(startIndex:endIndex) = 1.0e6_pReal ! Don't use for convergence check ! global alias plasticState(p)%slipRate => plasticState(p)%dotState(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 !-------------------------------------------------------------------------------------------------- -!> @brief calculates plastic velocity gradient and its tangent +!> @brief Calculate plastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- pure module subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp, & Mp,T,instance,of) @@ -312,7 +271,7 @@ pure module subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp, & Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp !< derivative of Lp with respect to the Mandel stress - + real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress real(pReal), intent(in) :: & @@ -320,18 +279,18 @@ pure module subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp, & 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) @@ -340,17 +299,17 @@ pure module subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp, & + 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 !-------------------------------------------------------------------------------------------------- -!> @brief calculates the rate of change of microstructure +!> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- module subroutine plastic_disloUCLA_dotState(Mp,T,instance,of) - + real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress real(pReal), intent(in) :: & @@ -358,7 +317,7 @@ module subroutine plastic_disloUCLA_dotState(Mp,T,instance,of) integer, intent(in) :: & instance, & of - + real(pReal) :: & VacancyDiffusion real(pReal), dimension(param(instance)%sum_N_sl) :: & @@ -369,16 +328,16 @@ module subroutine plastic_disloUCLA_dotState(Mp,T,instance,of) dot_rho_dip_formation, & dot_rho_dip_climb, & dip_distance - + 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 dot_rho_dip_formation = 0.0_pReal dot_rho_dip_climb = 0.0_pReal @@ -393,73 +352,73 @@ module subroutine plastic_disloUCLA_dotState(Mp,T,instance,of) * (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 - 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) = 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 - dot_rho_dip_climb - + end associate end subroutine plastic_disloUCLA_dotState !-------------------------------------------------------------------------------------------------- -!> @brief calculates derived quantities from state +!> @brief Calculate derived quantities from state. !-------------------------------------------------------------------------------------------------- module subroutine plastic_disloUCLA_dependentState(instance,of) - + integer, intent(in) :: & instance, & of - + real(pReal), dimension(param(instance)%sum_N_sl) :: & dislocationSpacing - + associate(prm => param(instance), stt => state(instance),dst => dependentState(instance)) - + dislocationSpacing = sqrt(matmul(prm%forestProjectionEdge, & stt%rho_mob(:,of)+stt%rho_dip(:,of))) dst%threshold_stress(:,of) = prm%mu*prm%b_sl & * sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,of)+stt%rho_dip(:,of))) - + dst%Lambda_sl(:,of) = prm%D/(1.0_pReal+prm%D*dislocationSpacing/prm%i_sl) - + end associate end subroutine plastic_disloUCLA_dependentState !-------------------------------------------------------------------------------------------------- -!> @brief writes results to HDF5 output file +!> @brief Write results to HDF5 output file. !-------------------------------------------------------------------------------------------------- module subroutine plastic_disloUCLA_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group - + integer :: o 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) - call results_writeDataset(group,stt%rho_mob,'rho_mob',& - 'mobile dislocation density','1/m²') - case (rho_dip_ID) - call results_writeDataset(group,stt%rho_dip,'rho_dip',& - 'dislocation dipole density''1/m²') - case (dot_gamma_sl_ID) - call results_writeDataset(group,stt%gamma_sl,'dot_gamma_sl',& ! this is not dot!! - 'plastic shear','1') - case (Lambda_sl_ID) - call results_writeDataset(group,dst%Lambda_sl,'Lambda_sl',& - 'mean free path for slip','m') - case (tau_pass_ID) - call results_writeDataset(group,dst%threshold_stress,'tau_pass',& - 'threshold stress for slip','Pa') + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) + case('edge_density') ! ToDo: should be rho_mob + if(prm%sum_N_sl>0) call results_writeDataset(group,stt%rho_mob,'rho_mob',& + 'mobile dislocation density','1/m²') + case('dipole_density') ! ToDo: should be rho_dip + if(prm%sum_N_sl>0) call results_writeDataset(group,stt%rho_dip,'rho_dip',& + 'dislocation dipole density''1/m²') + case('shear_rate_slip') ! should be gamma + if(prm%sum_N_sl>0) call results_writeDataset(group,stt%gamma_sl,'dot_gamma_sl',& ! this is not dot!! + 'plastic shear','1') + case('mfp_slip') !ToDo: should be Lambda + if(prm%sum_N_sl>0) call results_writeDataset(group,dst%Lambda_sl,'Lambda_sl',& + 'mean free path for slip','m') + case('threshold_stress_slip') !ToDo: should be tau_pass + if(prm%sum_N_sl>0) call results_writeDataset(group,dst%threshold_stress,'tau_pass',& + 'threshold stress for slip','Pa') end select enddo outputsLoop end associate @@ -468,15 +427,15 @@ end subroutine plastic_disloUCLA_results !-------------------------------------------------------------------------------------------------- -!> @brief Shear rates on slip systems, their derivatives with respect to resolved stress and the -! resolved stresss +!> @brief Calculate shear rates on slip systems, their derivatives with respect to resolved +! stress, and the resolved stress. !> @details Derivatives and resolved stress are calculated only optionally. ! 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,T,instance,of, & dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg,tau_pos_out,tau_neg_out) - + real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress real(pReal), intent(in) :: & @@ -484,7 +443,7 @@ pure subroutine kinetics(Mp,T,instance,of, & integer, intent(in) :: & instance, & of - + real(pReal), intent(out), dimension(param(instance)%sum_N_sl) :: & dot_gamma_pos, & dot_gamma_neg @@ -501,82 +460,82 @@ pure subroutine kinetics(Mp,T,instance,of, & 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) - + 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 = -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 = -1.0_pReal * t_k / tau_pos - + 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 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) - + 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 = -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 = -1.0_pReal * t_k / tau_neg - + 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 ddot_gamma_dtau_neg = 0.0_pReal end where significantNegativeTau2 end if - + end associate end associate diff --git a/src/constitutive_plastic_dislotwin.f90 b/src/constitutive_plastic_dislotwin.f90 index 94e0177a0..152c630ae 100644 --- a/src/constitutive_plastic_dislotwin.f90 +++ b/src/constitutive_plastic_dislotwin.f90 @@ -8,35 +8,17 @@ !> @details to be done !-------------------------------------------------------------------------------------------------- submodule(constitutive) plastic_dislotwin - - real(pReal), parameter :: & + + real(pReal), parameter :: & 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, & - resolved_stress_slip_ID, & - tau_pass_ID, & - edge_dipole_distance_ID, & - f_tw_ID, & - Lambda_tw_ID, & - resolved_stress_twin_ID, & - tau_hat_tw_ID, & - f_tr_ID - end enum - + type :: tParameters real(pReal) :: & mu, & nu, & D0, & !< prefactor for self-diffusion coefficient Qsd, & !< activation energy for dislocation climb - omega, & !< frequency factor for dislocation climb + omega, & !< frequency factor for dislocation climb D, & !< grain size p_sb, & !< p-exponent in shear band velocity q_sb, & !< q-exponent in shear band velocity @@ -58,8 +40,8 @@ submodule(constitutive) plastic_dislotwin aTol_f_tr, & !< absolute tolerance for integration of trans volume fraction gamma_fcc_hex, & !< Free energy difference between austensite and martensite i_tr, & !< - h !< Stack height of hex nucleus - real(pReal), dimension(:), allocatable :: & + h !< Stack height of hex nucleus + real(pReal), allocatable, dimension(:) :: & rho_mob_0, & !< initial unipolar dislocation density per slip system rho_dip_0, & !< initial dipole dislocation density per slip system b_sl, & !< absolute length of burgers vector [m] for each slip system @@ -79,40 +61,39 @@ submodule(constitutive) plastic_dislotwin s, & !< s-exponent in trans nucleation rate gamma_char, & !< characteristic shear for twins B !< drag coefficient - real(pReal), dimension(:,:), allocatable :: & - h_sl_sl, & !< + real(pReal), allocatable, dimension(:,:) :: & + h_sl_sl, & !< h_sl_tw, & !< - h_tw_tw, & !< - h_sl_tr, & !< - h_tr_tr !< - integer, dimension(:,:), allocatable :: & - fcc_twinNucleationSlipPair ! ToDo: Better name? Is also use for trans - real(pReal), dimension(:,:), allocatable :: & + h_tw_tw, & !< + h_sl_tr, & !< + h_tr_tr, & !< n0_sl, & !< slip system normal forestProjection, & C66 - real(pReal), dimension(:,:,:), allocatable :: & + real(pReal), allocatable, dimension(:,:,:) :: & P_tr, & P_sl, & P_tw, & C66_tw, & C66_tr - integer :: & + integer :: & sum_N_sl, & !< total number of active slip system sum_N_tw, & !< total number of active twin system - sum_N_tr !< total number of active transformation system - integer, dimension(:), allocatable :: & + sum_N_tr !< total number of active transformation system + integer, allocatable, dimension(:) :: & N_sl, & !< number of active slip systems for each family N_tw, & !< number of active twin systems for each family N_tr !< number of active transformation systems for each family - integer(kind(undefined_ID)), dimension(:), allocatable :: & - outputID !< ID of each post result output + integer, allocatable, dimension(:,:) :: & + fcc_twinNucleationSlipPair ! ToDo: Better name? Is also use for trans + character(len=pStringLen), allocatable, dimension(:) :: & + output logical :: & ExtendedDislocations, & !< consider split into partials for climb calculation fccTwinTransNucleation, & !< twinning and transformation models are for fcc dipoleFormation !< flag indicating consideration of dipole formation end type !< container type for internal constitutive parameters - + type :: tDislotwinState real(pReal), dimension(:,:), pointer :: & rho_mob, & @@ -121,7 +102,7 @@ submodule(constitutive) plastic_dislotwin f_tw, & f_tr end type tDislotwinState - + type :: tDislotwinMicrostructure real(pReal), dimension(:,:), allocatable :: & Lambda_sl, & !< mean free path between 2 obstacles seen by a moving dislocation @@ -135,7 +116,7 @@ submodule(constitutive) plastic_dislotwin tau_r_tw, & !< stress to bring partials close together (twin) tau_r_tr !< stress to bring partials close together (trans) end type tDislotwinMicrostructure - + !-------------------------------------------------------------------------------------------------- ! containers for parameters and state type(tParameters), allocatable, dimension(:) :: param @@ -148,7 +129,7 @@ contains !-------------------------------------------------------------------------------------------------- -!> @brief module initialization +!> @brief Perform module initialization. !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- module subroutine plastic_dislotwin_init @@ -156,39 +137,34 @@ module subroutine plastic_dislotwin_init integer :: & Ninstance, & p, i, & - NipcMyPhase, outputSize, & + NipcMyPhase, & sizeState, sizeDotState, & startIndex, endIndex - - integer(kind(undefined_ID)) :: & - outputID - + character(len=pStringLen) :: & extmsg = '' - character(len=pStringLen), dimension(:), allocatable :: & - outputs - - write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_DISLOTWIN_label//' init -+>>>' - - write(6,'(/,a)') ' Ma and Roters, Acta Materialia 52(12):3603–3612, 2004' - write(6,'(a)') ' https://doi.org/10.1016/j.actamat.2004.04.012' - - write(6,'(/,a)') ' Roters et al., Computational Materials Science 39:91–95, 2007' - write(6,'(a)') ' https://doi.org/10.1016/j.commatsci.2006.04.014' - - write(6,'(/,a)') ' Wong et al., Acta Materialia 118:140–151, 2016' - write(6,'(a,/)') ' https://doi.org/10.1016/j.actamat.2016.07.032' - + + write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_DISLOTWIN_label//' init -+>>>'; flush(6) + + write(6,'(/,a)') ' Ma and Roters, Acta Materialia 52(12):3603–3612, 2004' + write(6,'(a)') ' https://doi.org/10.1016/j.actamat.2004.04.012' + + write(6,'(/,a)') ' Roters et al., Computational Materials Science 39:91–95, 2007' + write(6,'(a)') ' https://doi.org/10.1016/j.commatsci.2006.04.014' + + write(6,'(/,a)') ' Wong et al., Acta Materialia 118:140–151, 2016' + write(6,'(a,/)') ' https://doi.org/10.1016/j.actamat.2016.07.032' + Ninstance = count(phase_plasticity == PLASTICITY_DISLOTWIN_ID) - + if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - + allocate(param(Ninstance)) allocate(state(Ninstance)) allocate(dotState(Ninstance)) allocate(dependentState(Ninstance)) - + do p = 1, size(phase_plasticity) if (phase_plasticity(p) /= PLASTICITY_DISLOTWIN_ID) cycle associate(prm => param(phase_plasticityInstance(p)), & @@ -196,17 +172,16 @@ module subroutine plastic_dislotwin_init stt => state(phase_plasticityInstance(p)), & dst => dependentState(phase_plasticityInstance(p)), & config => config_phase(p)) - + prm%aTol_rho = config%getFloat('atol_rho', defaultVal=0.0_pReal) prm%aTol_f_tw = config%getFloat('atol_twinfrac', defaultVal=0.0_pReal) prm%aTol_f_tr = config%getFloat('atol_transfrac', defaultVal=0.0_pReal) - + ! This data is read in already in lattice prm%mu = lattice_mu(p) prm%nu = lattice_nu(p) prm%C66 = lattice_C66(1:6,1:6,p) - - + !-------------------------------------------------------------------------------------------------- ! slip related parameters prm%N_sl = config%getInts('nslip',defaultVal=emptyIntArray) @@ -220,14 +195,14 @@ module subroutine plastic_dislotwin_init prm%forestProjection = lattice_forestProjection_edge(prm%N_sl,config%getString('lattice_structure'),& config%getFloat('c/a',defaultVal=0.0_pReal)) prm%forestProjection = transpose(prm%forestProjection) - + prm%n0_sl = lattice_slip_normal(prm%N_sl,config%getString('lattice_structure'),& - config%getFloat('c/a',defaultVal=0.0_pReal)) + config%getFloat('c/a',defaultVal=0.0_pReal)) prm%fccTwinTransNucleation = merge(.true., .false., lattice_structure(p) == LATTICE_FCC_ID) & .and. (prm%N_sl(1) == 12) if(prm%fccTwinTransNucleation) & prm%fcc_twinNucleationSlipPair = lattice_fcc_twinNucleationSlipPair - + 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)) @@ -238,7 +213,7 @@ module subroutine plastic_dislotwin_init prm%q = config%getFloats('q_slip', requiredSize=size(prm%N_sl)) prm%B = config%getFloats('b', requiredSize=size(prm%N_sl), & defaultVal=[(0.0_pReal, i=1,size(prm%N_sl))]) - + prm%tau_0 = config%getFloat('solidsolutionstrength') prm%CEdgeDipMinDistance = config%getFloat('cedgedipmindistance') prm%D0 = config%getFloat('d0') @@ -249,15 +224,15 @@ module subroutine plastic_dislotwin_init prm%SFE_0K = config%getFloat('sfe_0k') prm%dSFE_dT = config%getFloat('dsfe_dt') endif - + ! multiplication factor according to crystal structure (nearest neighbors bcc vs fcc/hex) !@details: Refer: Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981 prm%omega = config%getFloat('omega', defaultVal = 1000.0_pReal) & * merge(12.0_pReal, & 8.0_pReal, & lattice_structure(p) == LATTICE_FCC_ID .or. lattice_structure(p) == LATTICE_HEX_ID) - - + + ! 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) @@ -268,8 +243,8 @@ module subroutine plastic_dislotwin_init prm%p = math_expand(prm%p, prm%N_sl) prm%q = math_expand(prm%q, prm%N_sl) prm%B = math_expand(prm%B, prm%N_sl) - prm%atomicVolume = math_expand(prm%atomicVolume,prm%N_sl) - + prm%atomicVolume = math_expand(prm%atomicVolume,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' @@ -282,11 +257,11 @@ module subroutine plastic_dislotwin_init if (any(prm%B < 0.0_pReal)) extmsg = trim(extmsg)//' B' if (any(prm%p<=0.0_pReal .or. prm%p>1.0_pReal)) extmsg = trim(extmsg)//' p' if (any(prm%q< 1.0_pReal .or. prm%q>2.0_pReal)) extmsg = trim(extmsg)//' q' - + else slipActive allocate(prm%b_sl(0)) endif slipActive - + !-------------------------------------------------------------------------------------------------- ! twin related parameters prm%N_tw = config%getInts('ntwin', defaultVal=emptyIntArray) @@ -297,31 +272,31 @@ module subroutine plastic_dislotwin_init prm%h_tw_tw = lattice_interaction_TwinByTwin(prm%N_tw,& config%getFloats('interaction_twintwin'), & config%getString('lattice_structure')) - + prm%b_tw = config%getFloats('twinburgers', requiredSize=size(prm%N_tw)) prm%t_tw = config%getFloats('twinsize', requiredSize=size(prm%N_tw)) prm%r = config%getFloats('r_twin', requiredSize=size(prm%N_tw)) - + prm%xc_twin = config%getFloat('xc_twin') prm%L_tw = config%getFloat('l0_twin') prm%i_tw = config%getFloat('cmfptwin') - + prm%gamma_char= lattice_characteristicShear_Twin(prm%N_tw,config%getString('lattice_structure'),& config%getFloat('c/a',defaultVal=0.0_pReal)) - + prm%C66_tw = lattice_C66_twin(prm%N_tw,prm%C66,config%getString('lattice_structure'),& config%getFloat('c/a',defaultVal=0.0_pReal)) - + if (.not. prm%fccTwinTransNucleation) then - prm%dot_N_0_tw = config%getFloats('ndot0_twin') + prm%dot_N_0_tw = config%getFloats('ndot0_twin') prm%dot_N_0_tw = math_expand(prm%dot_N_0_tw,prm%N_tw) endif - + ! expand: family => system prm%b_tw = math_expand(prm%b_tw,prm%N_tw) prm%t_tw = math_expand(prm%t_tw,prm%N_tw) prm%r = math_expand(prm%r,prm%N_tw) - + else allocate(prm%gamma_char(0)) allocate(prm%t_tw (0)) @@ -329,7 +304,7 @@ module subroutine plastic_dislotwin_init allocate(prm%r (0)) allocate(prm%h_tw_tw (0,0)) endif - + !-------------------------------------------------------------------------------------------------- ! transformation related parameters prm%N_tr = config%getInts('ntrans', defaultVal=emptyIntArray) @@ -337,29 +312,29 @@ module subroutine plastic_dislotwin_init if (prm%sum_N_tr > 0) then prm%b_tr = config%getFloats('transburgers') prm%b_tr = math_expand(prm%b_tr,prm%N_tr) - + prm%h = config%getFloat('transstackheight', defaultVal=0.0_pReal) ! ToDo: How to handle that??? prm%i_tr = config%getFloat('cmfptrans', defaultVal=0.0_pReal) ! ToDo: How to handle that??? prm%gamma_fcc_hex = config%getFloat('deltag') prm%xc_trans = config%getFloat('xc_trans', defaultVal=0.0_pReal) ! ToDo: How to handle that??? prm%L_tr = config%getFloat('l0_trans') - + prm%h_tr_tr = lattice_interaction_TransByTrans(prm%N_tr,& config%getFloats('interaction_transtrans'), & config%getString('lattice_structure')) - + prm%C66_tr = lattice_C66_trans(prm%N_tr,prm%C66, & config%getString('trans_lattice_structure'), & 0.0_pReal, & config%getFloat('a_bcc', defaultVal=0.0_pReal), & config%getFloat('a_fcc', defaultVal=0.0_pReal)) - + prm%P_tr = lattice_SchmidMatrix_trans(prm%N_tr, & config%getString('trans_lattice_structure'), & 0.0_pReal, & config%getFloat('a_bcc', defaultVal=0.0_pReal), & config%getFloat('a_fcc', defaultVal=0.0_pReal)) - + if (lattice_structure(p) /= LATTICE_fcc_ID) then prm%dot_N_0_tr = config%getFloats('ndot0_trans') prm%dot_N_0_tr = math_expand(prm%dot_N_0_tr,prm%N_tr) @@ -374,59 +349,57 @@ module subroutine plastic_dislotwin_init allocate(prm%s (0)) allocate(prm%h_tr_tr(0,0)) endif - + if (sum(prm%N_tw) > 0 .or. prm%sum_N_tr > 0) then prm%SFE_0K = config%getFloat('sfe_0k') prm%dSFE_dT = config%getFloat('dsfe_dt') prm%V_cs = config%getFloat('vcrossslip') endif - + if (prm%sum_N_sl > 0 .and. prm%sum_N_tw > 0) then prm%h_sl_tw = lattice_interaction_SlipByTwin(prm%N_sl,prm%N_tw,& config%getFloats('interaction_sliptwin'), & config%getString('lattice_structure')) if (prm%fccTwinTransNucleation .and. prm%sum_N_tw > 12) write(6,*) 'mist' ! ToDo: implement better test. The model will fail also if N_tw is [6,6] - endif - - if (prm%sum_N_sl > 0 .and. prm%sum_N_tr > 0) then + endif + + if (prm%sum_N_sl > 0 .and. prm%sum_N_tr > 0) then prm%h_sl_tr = lattice_interaction_SlipByTrans(prm%N_sl,prm%N_tr,& config%getFloats('interaction_sliptrans'), & config%getString('lattice_structure')) if (prm%fccTwinTransNucleation .and. prm%sum_N_tr > 12) write(6,*) 'mist' ! ToDo: implement better test. The model will fail also if N_tr is [6,6] - endif - + endif + !-------------------------------------------------------------------------------------------------- ! shearband related parameters prm%sbVelocity = config%getFloat('shearbandvelocity',defaultVal=0.0_pReal) - if (prm%sbVelocity > 0.0_pReal) then + if (prm%sbVelocity > 0.0_pReal) then prm%sbResistance = config%getFloat('shearbandresistance') prm%sbQedge = config%getFloat('qedgepersbsystem') prm%p_sb = config%getFloat('p_shearband') prm%q_sb = config%getFloat('q_shearband') - + ! sanity checks if (prm%sbResistance < 0.0_pReal) extmsg = trim(extmsg)//' shearbandresistance' if (prm%sbQedge < 0.0_pReal) extmsg = trim(extmsg)//' qedgepersbsystem' if (prm%p_sb <= 0.0_pReal) extmsg = trim(extmsg)//' p_shearband' if (prm%q_sb <= 0.0_pReal) extmsg = trim(extmsg)//' q_shearband' endif - - - + prm%D = config%getFloat('grainsize') - + if (config%keyExists('dipoleformationfactor')) call IO_error(1,ext_msg='use /nodipoleformation/') prm%dipoleformation = .not. config%keyExists('/nodipoleformation/') - - + + !if (Ndot0PerTwinFamily(f,p) < 0.0_pReal) & ! call IO_error(211,el=p,ext_msg='dot_N_0_tw ('//PLASTICITY_DISLOTWIN_label//')') - + if (any(prm%atomicVolume <= 0.0_pReal)) & call IO_error(211,el=p,ext_msg='cAtomicVolume ('//PLASTICITY_DISLOTWIN_label//')') if (prm%sum_N_tw > 0) then if (prm%aTol_rho <= 0.0_pReal) & - call IO_error(211,el=p,ext_msg='aTol_rho ('//PLASTICITY_DISLOTWIN_label//')') + call IO_error(211,el=p,ext_msg='aTol_rho ('//PLASTICITY_DISLOTWIN_label//')') if (prm%aTol_f_tw <= 0.0_pReal) & call IO_error(211,el=p,ext_msg='aTol_f_tw ('//PLASTICITY_DISLOTWIN_label//')') endif @@ -434,50 +407,9 @@ module subroutine plastic_dislotwin_init if (prm%aTol_f_tr <= 0.0_pReal) & call IO_error(211,el=p,ext_msg='aTol_f_tr ('//PLASTICITY_DISLOTWIN_label//')') endif - - outputs = config%getStrings('(output)', defaultVal=emptyStringArray) - allocate(prm%outputID(0)) - do i= 1, size(outputs) - outputID = undefined_ID - select case(outputs(i)) - case ('rho_mob') - outputID = merge(rho_mob_ID,undefined_ID,prm%sum_N_sl > 0) - outputSize = prm%sum_N_sl - case ('rho_dip') - outputID = merge(rho_dip_ID,undefined_ID,prm%sum_N_sl > 0) - outputSize = prm%sum_N_sl - case ('gamma_sl') - outputID = merge(gamma_sl_ID,undefined_ID,prm%sum_N_sl > 0) - outputSize = prm%sum_N_sl - case ('lambda_sl') - outputID = merge(Lambda_sl_ID,undefined_ID,prm%sum_N_sl > 0) - outputSize = prm%sum_N_sl - case ('tau_pass') - outputID= merge(tau_pass_ID,undefined_ID,prm%sum_N_sl > 0) - outputSize = prm%sum_N_sl - - case ('f_tw') - outputID = merge(f_tw_ID,undefined_ID,prm%sum_N_tw >0) - outputSize = prm%sum_N_tw - case ('lambda_tw') - outputID = merge(Lambda_tw_ID,undefined_ID,prm%sum_N_tw >0) - outputSize = prm%sum_N_tw - case ('tau_hat_tw') - outputID = merge(tau_hat_tw_ID,undefined_ID,prm%sum_N_tw >0) - outputSize = prm%sum_N_tw - - case ('f_tr') - outputID = f_tr_ID - outputSize = prm%sum_N_tr - - end select - - if (outputID /= undefined_ID) then - prm%outputID = [prm%outputID, outputID] - endif - - enddo - + + prm%output = config%getStrings('(output)', defaultVal=emptyStringArray) + !-------------------------------------------------------------------------------------------------- ! allocate state arrays NipcMyPhase = count(material_phaseAt == p) * discretization_nIP @@ -485,9 +417,9 @@ module subroutine plastic_dislotwin_init + size(['f_tw']) * prm%sum_N_tw & + size(['f_tr']) * prm%sum_N_tr sizeState = sizeDotState - + call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0) - + !-------------------------------------------------------------------------------------------------- ! locally defined state aliases and initialization of state0 and aTolState startIndex = 1 @@ -496,14 +428,14 @@ module subroutine plastic_dislotwin_init 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,:) @@ -511,66 +443,66 @@ module subroutine plastic_dislotwin_init plasticState(p)%aTolState(startIndex:endIndex) = 1.0e6_pReal !ToDo: better make optional parameter ! global alias plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) - + startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_tw stt%f_tw=>plasticState(p)%state(startIndex:endIndex,:) dot%f_tw=>plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTol_f_tw - + startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_tr stt%f_tr=>plasticState(p)%state(startIndex:endIndex,:) dot%f_tr=>plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTol_f_tr - + allocate(dst%Lambda_sl (prm%sum_N_sl,NipcMyPhase),source=0.0_pReal) allocate(dst%tau_pass (prm%sum_N_sl,NipcMyPhase),source=0.0_pReal) - + allocate(dst%Lambda_tw (prm%sum_N_tw,NipcMyPhase),source=0.0_pReal) allocate(dst%tau_hat_tw (prm%sum_N_tw,NipcMyPhase),source=0.0_pReal) allocate(dst%tau_r_tw (prm%sum_N_tw,NipcMyPhase),source=0.0_pReal) allocate(dst%V_tw (prm%sum_N_tw,NipcMyPhase),source=0.0_pReal) - + allocate(dst%Lambda_tr (prm%sum_N_tr,NipcMyPhase),source=0.0_pReal) allocate(dst%tau_hat_tr (prm%sum_N_tr,NipcMyPhase),source=0.0_pReal) allocate(dst%tau_r_tr (prm%sum_N_tr,NipcMyPhase),source=0.0_pReal) allocate(dst%V_tr (prm%sum_N_tr,NipcMyPhase),source=0.0_pReal) - - + + plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally - + end associate - + enddo end subroutine plastic_dislotwin_init !-------------------------------------------------------------------------------------------------- -!> @brief returns the homogenized elasticity matrix +!> @brief Return the homogenized elasticity matrix. !-------------------------------------------------------------------------------------------------- module function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC) - + real(pReal), dimension(6,6) :: & homogenizedC integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element - + integer :: i, & of real(pReal) :: f_unrotated - + of = material_phasememberAt(ipc,ip,el) associate(prm => param(phase_plasticityInstance(material_phaseAt(ipc,el))),& stt => state(phase_plasticityInstance(material_phaseAT(ipc,el)))) - + f_unrotated = 1.0_pReal & - sum(stt%f_tw(1:prm%sum_N_tw,of)) & - sum(stt%f_tr(1:prm%sum_N_tr,of)) - + homogenizedC = f_unrotated * prm%C66 do i=1,prm%sum_N_tw homogenizedC = homogenizedC & @@ -580,23 +512,23 @@ module function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC) homogenizedC = homogenizedC & + stt%f_tr(i,of)*prm%C66_tr(1:6,1:6,i) enddo - + end associate - + end function plastic_dislotwin_homogenizedC !-------------------------------------------------------------------------------------------------- -!> @brief calculates plastic velocity gradient and its tangent +!> @brief Calculate plastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) - + real(pReal), dimension(3,3), intent(out) :: Lp real(pReal), dimension(3,3,3,3), intent(out) :: dLp_dMp real(pReal), dimension(3,3), intent(in) :: Mp integer, intent(in) :: instance,of real(pReal), intent(in) :: T - + integer :: i,k,l,m,n real(pReal) :: & f_unrotated,StressRatio_p,& @@ -632,16 +564,16 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) 0, 1,-1, & 0, 1, 1 & ],pReal),[ 3,6]) - + associate(prm => param(instance), stt => state(instance)) - + f_unrotated = 1.0_pReal & - sum(stt%f_tw(1:prm%sum_N_tw,of)) & - sum(stt%f_tr(1:prm%sum_N_tr,of)) - + Lp = 0.0_pReal - dLp_dMp = 0.0_pReal - + dLp_dMp = 0.0_pReal + call kinetics_slip(Mp,T,instance,of,dot_gamma_sl,ddot_gamma_dtau_slip) slipContribution: do i = 1, prm%sum_N_sl Lp = Lp + dot_gamma_sl(i)*prm%P_sl(1:3,1:3,i) @@ -649,37 +581,37 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + ddot_gamma_dtau_slip(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i) enddo slipContribution - + !ToDo: Why do this before shear banding? Lp = Lp * f_unrotated dLp_dMp = dLp_dMp * f_unrotated - + shearBandingContribution: if(dNeq0(prm%sbVelocity)) then - + BoltzmannRatio = prm%sbQedge/(kB*T) call math_eigenValuesVectorsSym(Mp,eigValues,eigVectors,error) - + do i = 1,6 P_sb = 0.5_pReal * math_outer(matmul(eigVectors,sb_sComposition(1:3,i)),& matmul(eigVectors,sb_mComposition(1:3,i))) tau = math_mul33xx33(Mp,P_sb) - + significantShearBandStress: if (abs(tau) > tol_math_check) then StressRatio_p = (abs(tau)/prm%sbResistance)**prm%p_sb dot_gamma_sb = sign(prm%sbVelocity*exp(-BoltzmannRatio*(1-StressRatio_p)**prm%q_sb), tau) ddot_gamma_dtau = abs(dot_gamma_sb)*BoltzmannRatio* prm%p_sb*prm%q_sb/ prm%sbResistance & * (abs(tau)/prm%sbResistance)**(prm%p_sb-1.0_pReal) & * (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal) - + Lp = Lp + dot_gamma_sb * P_sb 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 * P_sb(k,l) * P_sb(m,n) endif significantShearBandStress enddo - + endif shearBandingContribution - + call kinetics_twin(Mp,T,dot_gamma_sl,instance,of,dot_gamma_twin,ddot_gamma_dtau_twin) twinContibution: do i = 1, prm%sum_N_tw Lp = Lp + dot_gamma_twin(i)*prm%P_tw(1:3,1:3,i) * f_unrotated @@ -687,7 +619,7 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + ddot_gamma_dtau_twin(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i) * f_unrotated enddo twinContibution - + call kinetics_trans(Mp,T,dot_gamma_sl,instance,of,dot_gamma_tr,ddot_gamma_dtau_trans) transContibution: do i = 1, prm%sum_N_tr Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i) * f_unrotated @@ -695,15 +627,15 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + ddot_gamma_dtau_trans(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i) * f_unrotated enddo transContibution - - + + end associate - + end subroutine plastic_dislotwin_LpAndItsTangent !-------------------------------------------------------------------------------------------------- -!> @brief calculates the rate of change of microstructure +!> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) @@ -720,10 +652,10 @@ module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) f_unrotated, & VacancyDiffusion, & rho_dip_distance, & - v_cl, & !< climb velocity + v_cl, & !< climb velocity Gamma, & !< stacking fault energy tau, & - sigma_cl, & !< climb stress + sigma_cl, & !< climb stress b_d !< ratio of burgers vector to stacking fault width real(pReal), dimension(param(instance)%sum_N_sl) :: & dot_rho_dip_formation, & @@ -745,9 +677,9 @@ module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) call kinetics_slip(Mp,T,instance,of,dot_gamma_sl) dot%gamma_sl(:,of) = abs(dot_gamma_sl) - + rho_dip_distance_min = prm%CEdgeDipMinDistance*prm%b_sl - + slipState: do i = 1, prm%sum_N_sl tau = math_mul33xx33(Mp,prm%P_sl(1:3,1:3,i)) @@ -771,15 +703,15 @@ module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) else !@details: Refer: Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981 sigma_cl = dot_product(prm%n0_sl(1:3,i),matmul(Mp,prm%n0_sl(1:3,i))) - if (prm%ExtendedDislocations) then + if (prm%ExtendedDislocations) then Gamma = prm%SFE_0K + prm%dSFE_dT * T b_d = 24.0_pReal*PI*(1.0_pReal - prm%nu)/(2.0_pReal + prm%nu)* Gamma/(prm%mu*prm%b_sl(i)) else - b_d = 1.0_pReal + b_d = 1.0_pReal endif v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Qsd/(kB*T)) & * (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(kB*T)) - 1.0_pReal) - + dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,of) & / (rho_dip_distance-rho_dip_distance_min(i)) endif @@ -801,12 +733,12 @@ module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) dot%f_tr(:,of) = f_unrotated*dot_gamma_tr end associate - + end subroutine plastic_dislotwin_dotState !-------------------------------------------------------------------------------------------------- -!> @brief calculates derived quantities from state +!> @brief Calculate derived quantities from state. !-------------------------------------------------------------------------------------------------- module subroutine plastic_dislotwin_dependentState(T,instance,of) @@ -815,7 +747,7 @@ module subroutine plastic_dislotwin_dependentState(T,instance,of) of real(pReal), intent(in) :: & T - + real(pReal) :: & sumf_twin,Gamma,sumf_trans real(pReal), dimension(param(instance)%sum_N_sl) :: & @@ -830,35 +762,35 @@ module subroutine plastic_dislotwin_dependentState(T,instance,of) f_over_t_tr real(pReal), dimension(:), allocatable :: & x0 - - + + associate(prm => param(instance),& stt => state(instance),& dst => dependentState(instance)) - + sumf_twin = sum(stt%f_tw(1:prm%sum_N_tw,of)) sumf_trans = sum(stt%f_tr(1:prm%sum_N_tr,of)) - + Gamma = prm%SFE_0K + prm%dSFE_dT * T - + !* rescaled volume fraction for topology f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,of)/prm%t_tw ! this is per system ... f_over_t_tr = sumf_trans/prm%t_tr ! but this not ! ToDo ...Physically correct, but naming could be adjusted - + inv_lambda_sl_sl = sqrt(matmul(prm%forestProjection, & stt%rho_mob(:,of)+stt%rho_dip(:,of)))/prm%CLambdaSlip - + if (prm%sum_N_tw > 0 .and. prm%sum_N_sl > 0) & inv_lambda_sl_tw = matmul(prm%h_sl_tw,f_over_t_tw)/(1.0_pReal-sumf_twin) - + inv_lambda_tw_tw = matmul(prm%h_tw_tw,f_over_t_tw)/(1.0_pReal-sumf_twin) - + if (prm%sum_N_tr > 0 .and. prm%sum_N_sl > 0) & inv_lambda_sl_tr = matmul(prm%h_sl_tr,f_over_t_tr)/(1.0_pReal-sumf_trans) - + inv_lambda_tr_tr = matmul(prm%h_tr_tr,f_over_t_tr)/(1.0_pReal-sumf_trans) - + if ((prm%sum_N_tw > 0) .or. (prm%sum_N_tr > 0)) then ! ToDo: better logic needed here dst%Lambda_sl(:,of) = prm%D & / (1.0_pReal+prm%D*(inv_lambda_sl_sl + inv_lambda_sl_tw + inv_lambda_sl_tr)) @@ -866,13 +798,13 @@ module subroutine plastic_dislotwin_dependentState(T,instance,of) dst%Lambda_sl(:,of) = prm%D & / (1.0_pReal+prm%D*inv_lambda_sl_sl) !!!!!! correct? endif - + dst%Lambda_tw(:,of) = prm%i_tw*prm%D/(1.0_pReal+prm%D*inv_lambda_tw_tw) dst%Lambda_tr(:,of) = prm%i_tr*prm%D/(1.0_pReal+prm%D*inv_lambda_tr_tr) - + !* threshold stress for dislocation motion dst%tau_pass(:,of) = prm%mu*prm%b_sl* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,of)+stt%rho_dip(:,of))) - + !* threshold stress for growing twin/martensite if(prm%sum_N_tw == prm%sum_N_sl) & dst%tau_hat_tw(:,of) = Gamma/(3.0_pReal*prm%b_tw) & @@ -880,66 +812,67 @@ module subroutine plastic_dislotwin_dependentState(T,instance,of) if(prm%sum_N_tr == prm%sum_N_sl) & dst%tau_hat_tr(:,of) = Gamma/(3.0_pReal*prm%b_tr) & + 3.0_pReal*prm%b_tr*prm%mu/(prm%L_tr*prm%b_sl) & ! slip burgers here correct? - + prm%h*prm%gamma_fcc_hex/ (3.0_pReal*prm%b_tr) - + + prm%h*prm%gamma_fcc_hex/ (3.0_pReal*prm%b_tr) + dst%V_tw(:,of) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,of)**2.0_pReal dst%V_tr(:,of) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,of)**2.0_pReal - - + + x0 = prm%mu*prm%b_tw**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the burgers vector for slip and is the same for twin and trans dst%tau_r_tw(:,of) = prm%mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%xc_twin)+cos(pi/3.0_pReal)/x0) - + x0 = prm%mu*prm%b_tr**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the burgers vector for slip dst%tau_r_tr(:,of) = prm%mu*prm%b_tr/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%xc_trans)+cos(pi/3.0_pReal)/x0) - + end associate end subroutine plastic_dislotwin_dependentState !-------------------------------------------------------------------------------------------------- -!> @brief writes results to HDF5 output file +!> @brief Write results to HDF5 output file. !-------------------------------------------------------------------------------------------------- module subroutine plastic_dislotwin_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group + integer :: o - associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) - outputsLoop: do o = 1,size(prm%outputID) - select case(prm%outputID(o)) + associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) - case (rho_mob_ID) - call results_writeDataset(group,stt%rho_mob,'rho_mob',& - 'mobile dislocation density','1/m²') - case (rho_dip_ID) - call results_writeDataset(group,stt%rho_dip,'rho_dip',& - 'dislocation dipole density''1/m²') - case (gamma_sl_ID) - call results_writeDataset(group,stt%gamma_sl,'gamma_sl',& - 'plastic shear','1') - case (Lambda_sl_ID) - call results_writeDataset(group,dst%Lambda_sl,'Lambda_sl',& - 'mean free path for slip','m') - case (tau_pass_ID) - call results_writeDataset(group,dst%tau_pass,'tau_pass',& - 'passing stress for slip','Pa') + case('rho_mob') + if(prm%sum_N_sl>0) call results_writeDataset(group,stt%rho_mob,'rho_mob',& + 'mobile dislocation density','1/m²') + case('rho_dip') + if(prm%sum_N_sl>0) call results_writeDataset(group,stt%rho_dip,'rho_dip',& + 'dislocation dipole density''1/m²') + case('gamma_sl') + if(prm%sum_N_sl>0) call results_writeDataset(group,stt%gamma_sl,'gamma_sl',& + 'plastic shear','1') + case('lambda_sl') + if(prm%sum_N_sl>0) call results_writeDataset(group,dst%Lambda_sl,'Lambda_sl',& + 'mean free path for slip','m') + case('tau_pass') + if(prm%sum_N_sl>0) call results_writeDataset(group,dst%tau_pass,'tau_pass',& + 'passing stress for slip','Pa') - case (f_tw_ID) - call results_writeDataset(group,stt%f_tw,'f_tw',& - 'twinned volume fraction','m³/m³') - case (Lambda_tw_ID) - call results_writeDataset(group,dst%Lambda_tw,'Lambda_tw',& - 'mean free path for twinning','m') - case (tau_hat_tw_ID) - call results_writeDataset(group,dst%tau_hat_tw,'tau_hat_tw',& - 'threshold stress for twinning','Pa') + case('f_tw') + if(prm%sum_N_tw>0) call results_writeDataset(group,stt%f_tw,'f_tw',& + 'twinned volume fraction','m³/m³') + case('lambda_tw') + if(prm%sum_N_tw>0) call results_writeDataset(group,dst%Lambda_tw,'Lambda_tw',& + 'mean free path for twinning','m') + case('tau_hat_tw') + if(prm%sum_N_tw>0) call results_writeDataset(group,dst%tau_hat_tw,'tau_hat_tw',& + 'threshold stress for twinning','Pa') + + case('f_tr') + if(prm%sum_N_tr>0) call results_writeDataset(group,stt%f_tr,'f_tr',& + 'martensite volume fraction','m³/m³') - case (f_tr_ID) - call results_writeDataset(group,stt%f_tr,'f_tr',& - 'martensite volume fraction','m³/m³') - end select enddo outputsLoop end associate @@ -948,8 +881,8 @@ end subroutine plastic_dislotwin_results !-------------------------------------------------------------------------------------------------- -!> @brief Shear rates on slip systems, their derivatives with respect to resolved stress and the -! resolved stresss +!> @brief Calculate shear rates on slip systems, their derivatives with respect to resolved +! stress, and the resolved stress. !> @details Derivatives and resolved stress are calculated only optionally. ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end @@ -964,7 +897,7 @@ pure subroutine kinetics_slip(Mp,T,instance,of, & integer, intent(in) :: & instance, & of - + real(pReal), dimension(param(instance)%sum_N_sl), intent(out) :: & dot_gamma_sl real(pReal), dimension(param(instance)%sum_N_sl), optional, intent(out) :: & @@ -972,7 +905,7 @@ pure subroutine kinetics_slip(Mp,T,instance,of, & tau_slip real(pReal), dimension(param(instance)%sum_N_sl) :: & ddot_gamma_dtau - + real(pReal), dimension(param(instance)%sum_N_sl) :: & tau, & stressRatio, & @@ -984,25 +917,25 @@ pure subroutine kinetics_slip(Mp,T,instance,of, & dV_run_inverse_dTau, & dV_dTau, & tau_eff !< effective resolved stress - integer :: i - + integer :: i + associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) - + do i = 1, prm%sum_N_sl tau(i) = math_mul33xx33(Mp,prm%P_sl(1:3,1:3,i)) enddo - + tau_eff = abs(tau)-dst%tau_pass(:,of) - + significantStress: where(tau_eff > tol_math_check) stressRatio = tau_eff/prm%tau_0 StressRatio_p = stressRatio** prm%p BoltzmannRatio = prm%Delta_F/(kB*T) v_wait_inverse = prm%v0**(-1.0_pReal) * exp(BoltzmannRatio*(1.0_pReal-StressRatio_p)** prm%q) v_run_inverse = prm%B/(tau_eff*prm%b_sl) - + dot_gamma_sl = sign(stt%rho_mob(:,of)*prm%b_sl/(v_wait_inverse+v_run_inverse),tau) - + dV_wait_inverse_dTau = -1.0_pReal * v_wait_inverse * prm%p * prm%q * BoltzmannRatio & * (stressRatio**(prm%p-1.0_pReal)) & * (1.0_pReal-StressRatio_p)**(prm%q-1.0_pReal) & @@ -1015,17 +948,21 @@ pure subroutine kinetics_slip(Mp,T,instance,of, & dot_gamma_sl = 0.0_pReal ddot_gamma_dtau = 0.0_pReal end where significantStress - + end associate - + if(present(ddot_gamma_dtau_slip)) ddot_gamma_dtau_slip = ddot_gamma_dtau if(present(tau_slip)) tau_slip = tau - + end subroutine kinetics_slip !-------------------------------------------------------------------------------------------------- -!> @brief calculates shear rates on twin systems +!> @brief Calculate shear rates on twin systems and their derivatives with respect to resolved +! stress. +!> @details Derivatives are calculated only optionally. +! 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_twin(Mp,T,dot_gamma_sl,instance,of,& dot_gamma_twin,ddot_gamma_dtau_twin) @@ -1039,22 +976,22 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,of,& of real(pReal), dimension(param(instance)%sum_N_sl), intent(in) :: & dot_gamma_sl - + real(pReal), dimension(param(instance)%sum_N_tw), intent(out) :: & dot_gamma_twin real(pReal), dimension(param(instance)%sum_N_tw), optional, intent(out) :: & ddot_gamma_dtau_twin - + real, dimension(param(instance)%sum_N_tw) :: & tau, & Ndot0, & stressRatio_r, & ddot_gamma_dtau - + integer :: i,s1,s2 - + associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) - + do i = 1, prm%sum_N_tw tau(i) = math_mul33xx33(Mp,prm%P_tw(1:3,1:3,i)) isFCC: if (prm%fccTwinTransNucleation) then @@ -1072,7 +1009,7 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,of,& Ndot0=prm%dot_N_0_tw(i) endif isFCC enddo - + significantStress: where(tau > tol_math_check) StressRatio_r = (dst%tau_hat_tw(:,of)/tau)**prm%r dot_gamma_twin = prm%gamma_char * dst%V_tw(:,of) * Ndot0*exp(-StressRatio_r) @@ -1081,16 +1018,20 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,of,& dot_gamma_twin = 0.0_pReal ddot_gamma_dtau = 0.0_pReal end where significantStress - + end associate - + if(present(ddot_gamma_dtau_twin)) ddot_gamma_dtau_twin = ddot_gamma_dtau end subroutine kinetics_twin !-------------------------------------------------------------------------------------------------- -!> @brief calculates shear rates on twin systems +!> @brief Calculate shear rates on transformation systems and their derivatives with respect to +! resolved stress. +!> @details Derivatives are calculated only optionally. +! 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_trans(Mp,T,dot_gamma_sl,instance,of,& dot_gamma_tr,ddot_gamma_dtau_trans) @@ -1104,21 +1045,21 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,of,& of real(pReal), dimension(param(instance)%sum_N_sl), intent(in) :: & dot_gamma_sl - + real(pReal), dimension(param(instance)%sum_N_tr), intent(out) :: & dot_gamma_tr real(pReal), dimension(param(instance)%sum_N_tr), optional, intent(out) :: & ddot_gamma_dtau_trans - + real, dimension(param(instance)%sum_N_tr) :: & tau, & Ndot0, & stressRatio_s, & ddot_gamma_dtau - + integer :: i,s1,s2 associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) - + do i = 1, prm%sum_N_tr tau(i) = math_mul33xx33(Mp,prm%P_tr(1:3,1:3,i)) isFCC: if (prm%fccTwinTransNucleation) then @@ -1136,7 +1077,7 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,of,& Ndot0=prm%dot_N_0_tr(i) endif isFCC enddo - + significantStress: where(tau > tol_math_check) StressRatio_s = (dst%tau_hat_tr(:,of)/tau)**prm%s dot_gamma_tr = dst%V_tr(:,of) * Ndot0*exp(-StressRatio_s) @@ -1145,9 +1086,9 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,of,& dot_gamma_tr = 0.0_pReal ddot_gamma_dtau = 0.0_pReal end where significantStress - + end associate - + if(present(ddot_gamma_dtau_trans)) ddot_gamma_dtau_trans = ddot_gamma_dtau end subroutine kinetics_trans diff --git a/src/constitutive_plastic_isotropic.f90 b/src/constitutive_plastic_isotropic.f90 index ef8a823d4..beed2112a 100644 --- a/src/constitutive_plastic_isotropic.f90 +++ b/src/constitutive_plastic_isotropic.f90 @@ -8,14 +8,7 @@ !! untextured polycrystal !-------------------------------------------------------------------------------------------------- submodule(constitutive) plastic_isotropic - - enum, bind(c) - enumerator :: & - undefined_ID, & - xi_ID, & - dot_gamma_ID - end enum - + type :: tParameters real(pReal) :: & M, & !< Taylor factor @@ -34,12 +27,12 @@ submodule(constitutive) plastic_isotropic aTol_gamma integer :: & of_debug = 0 - integer(kind(undefined_ID)), allocatable, dimension(:) :: & - outputID logical :: & dilatation + character(len=pStringLen), allocatable, dimension(:) :: & + output end type tParameters - + type :: tIsotropicState real(pReal), pointer, dimension(:) :: & xi, & @@ -56,50 +49,45 @@ submodule(constitutive) plastic_isotropic contains !-------------------------------------------------------------------------------------------------- -!> @brief module initialization +!> @brief Perform module initialization. !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- module subroutine plastic_isotropic_init integer :: & Ninstance, & - p, i, & + p, & NipcMyPhase, & sizeState, sizeDotState - - integer(kind(undefined_ID)) :: & - outputID - + character(len=pStringLen) :: & extmsg = '' - character(len=pStringLen), dimension(:), allocatable :: & - outputs - - write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_ISOTROPIC_label//' init -+>>>' - - write(6,'(/,a)') ' Maiti and Eisenlohr, Scripta Materialia 145:37–40, 2018' - write(6,'(a)') ' https://doi.org/10.1016/j.scriptamat.2017.09.047' - + + write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_ISOTROPIC_label//' init -+>>>'; flush(6) + + write(6,'(/,a)') ' Maiti and Eisenlohr, Scripta Materialia 145:37–40, 2018' + write(6,'(a)') ' https://doi.org/10.1016/j.scriptamat.2017.09.047' + Ninstance = count(phase_plasticity == PLASTICITY_ISOTROPIC_ID) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - + allocate(param(Ninstance)) allocate(state(Ninstance)) allocate(dotState(Ninstance)) - + do p = 1, size(phase_plasticity) if (phase_plasticity(p) /= PLASTICITY_ISOTROPIC_ID) cycle associate(prm => param(phase_plasticityInstance(p)), & dot => dotState(phase_plasticityInstance(p)), & stt => state(phase_plasticityInstance(p)), & config => config_phase(p)) - + #ifdef DEBUG if (p==material_phaseAt(debug_g,debug_e)) & prm%of_debug = material_phasememberAt(debug_g,debug_i,debug_e) #endif - + prm%xi_0 = config%getFloat('tau0') prm%xi_inf = config%getFloat('tausat') prm%dot_gamma_0 = config%getFloat('gdot0') @@ -114,9 +102,9 @@ module subroutine plastic_isotropic_init prm%a = config%getFloat('a') prm%aTol_xi = config%getFloat('atol_flowstress',defaultVal=1.0_pReal) prm%aTol_gamma = config%getFloat('atol_shear', defaultVal=1.0e-6_pReal) - + prm%dilatation = config%keyExists('/dilatation/') - + !-------------------------------------------------------------------------------------------------- ! sanity checks extmsg = '' @@ -128,79 +116,62 @@ module subroutine plastic_isotropic_init if (prm%M <= 0.0_pReal) extmsg = trim(extmsg)//' m' if (prm%aTol_xi <= 0.0_pReal) extmsg = trim(extmsg)//' atol_xi' if (prm%aTol_gamma <= 0.0_pReal) extmsg = trim(extmsg)//' atol_shear' - + !-------------------------------------------------------------------------------------------------- ! exit if any parameter is out of range if (extmsg /= '') & call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_ISOTROPIC_label//')') - + !-------------------------------------------------------------------------------------------------- ! output pararameters - outputs = config%getStrings('(output)',defaultVal=emptyStringArray) - allocate(prm%outputID(0)) - do i=1, size(outputs) - outputID = undefined_ID - select case(outputs(i)) - - case ('flowstress') - outputID = xi_ID - case ('strainrate') - outputID = dot_gamma_ID - - end select - - if (outputID /= undefined_ID) then - prm%outputID = [prm%outputID, outputID] - endif - - enddo - + prm%output = config%getStrings('(output)',defaultVal=emptyStringArray) + !-------------------------------------------------------------------------------------------------- ! allocate state arrays NipcMyPhase = count(material_phaseAt == p) * discretization_nIP sizeDotState = size(['xi ','accumulated_shear']) sizeState = sizeDotState - + call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0) - + !-------------------------------------------------------------------------------------------------- ! locally defined state aliases and initialization of state0 and aTolState stt%xi => plasticState(p)%state (1,:) stt%xi = prm%xi_0 dot%xi => plasticState(p)%dotState(1,:) plasticState(p)%aTolState(1) = prm%aTol_xi - + stt%gamma => plasticState(p)%state (2,:) dot%gamma => plasticState(p)%dotState(2,:) plasticState(p)%aTolState(2) = prm%aTol_gamma ! global alias plasticState(p)%slipRate => plasticState(p)%dotState(2:2,:) - + plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally - + end associate - + enddo end subroutine plastic_isotropic_init !-------------------------------------------------------------------------------------------------- -!> @brief calculates plastic velocity gradient and its tangent +!> @brief Calculate plastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) - + real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp !< derivative of Lp with respect to the Mandel stress - + real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & instance, & of - + real(pReal), dimension(3,3) :: & Mp_dev !< deviatoric part of the Mandel stress real(pReal) :: & @@ -209,16 +180,16 @@ module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) squarenorm_Mp_dev !< square of the norm of the deviatoric part of the Mandel stress integer :: & k, l, m, n - + associate(prm => param(instance), stt => state(instance)) - + Mp_dev = math_deviatoric33(Mp) squarenorm_Mp_dev = math_mul33xx33(Mp_dev,Mp_dev) norm_Mp_dev = sqrt(squarenorm_Mp_dev) - + if (norm_Mp_dev > 0.0_pReal) then dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%M*stt%xi(of))) **prm%n - + Lp = dot_gamma/prm%M * Mp_dev/norm_Mp_dev #ifdef DEBUG if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0 & @@ -240,42 +211,42 @@ module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) Lp = 0.0_pReal dLp_dMp = 0.0_pReal end if - + end associate end subroutine plastic_isotropic_LpAndItsTangent !-------------------------------------------------------------------------------------------------- -!> @brief calculates plastic velocity gradient and its tangent +!> @brief Calculate inelastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of) - + real(pReal), dimension(3,3), intent(out) :: & Li !< inleastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLi_dMi !< derivative of Li with respect to Mandel stress - + real(pReal), dimension(3,3), intent(in) :: & - Mi !< Mandel stress + Mi !< Mandel stress integer, intent(in) :: & instance, & of - + real(pReal) :: & tr !< trace of spherical part of Mandel stress (= 3 x pressure) integer :: & k, l, m, n - + associate(prm => param(instance), stt => state(instance)) - + tr=math_trace33(math_spherical33(Mi)) if (prm%dilatation .and. abs(tr) > 0.0_pReal) then ! no stress or J2 plasticity --> Li and its derivative are zero Li = math_I3 & * prm%dot_gamma_0/prm%M * (3.0_pReal*prm%M*stt%xi(of))**(-prm%n) & * tr * abs(tr)**(prm%n-1.0_pReal) - + #ifdef DEBUG if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0 & .and. (of == prm%of_debug .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0)) then @@ -292,38 +263,38 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of) Li = 0.0_pReal dLi_dMi = 0.0_pReal endif - + end associate end subroutine plastic_isotropic_LiAndItsTangent !-------------------------------------------------------------------------------------------------- -!> @brief calculates the rate of change of microstructure +!> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- module subroutine plastic_isotropic_dotState(Mp,instance,of) - + real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & instance, & of - + real(pReal) :: & dot_gamma, & !< strainrate xi_inf_star, & !< saturation xi norm_Mp !< norm of the (deviatoric) Mandel stress - + associate(prm => param(instance), stt => state(instance), dot => dotState(instance)) - + if (prm%dilatation) then norm_Mp = sqrt(math_mul33xx33(Mp,Mp)) else norm_Mp = sqrt(math_mul33xx33(math_deviatoric33(Mp),math_deviatoric33(Mp))) endif - + dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp /(prm%M*stt%xi(of))) **prm%n - + if (dot_gamma > 1e-12_pReal) then if (dEq0(prm%c_1)) then xi_inf_star = prm%xi_inf @@ -339,28 +310,28 @@ module subroutine plastic_isotropic_dotState(Mp,instance,of) else dot%xi(of) = 0.0_pReal endif - + dot%gamma(of) = dot_gamma ! ToDo: not really used - + end associate end subroutine plastic_isotropic_dotState !-------------------------------------------------------------------------------------------------- -!> @brief writes results to HDF5 output file +!> @brief Write results to HDF5 output file. !-------------------------------------------------------------------------------------------------- module subroutine plastic_isotropic_results(instance,group) - integer, intent(in) :: instance + integer, intent(in) :: instance character(len=*), intent(in) :: group - + integer :: o associate(prm => param(instance), stt => state(instance)) - outputsLoop: do o = 1,size(prm%outputID) - select case(prm%outputID(o)) - case (xi_ID) + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) + case ('flowstress') ! ToDo: should be 'xi' call results_writeDataset(group,stt%xi,'xi','resistance against plastic flow','Pa') end select enddo outputsLoop diff --git a/src/constitutive_plastic_kinehardening.f90 b/src/constitutive_plastic_kinehardening.f90 index 11ce7c494..65751989c 100644 --- a/src/constitutive_plastic_kinehardening.f90 +++ b/src/constitutive_plastic_kinehardening.f90 @@ -7,26 +7,13 @@ !-------------------------------------------------------------------------------------------------- submodule(constitutive) plastic_kinehardening - enum, bind(c) - enumerator :: & - undefined_ID, & - crss_ID, & !< critical resolved stress - crss_back_ID, & !< critical resolved back stress - sense_ID, & !< sense of acting shear stress (-1 or +1) - chi0_ID, & !< backstress at last switch of stress sense (positive?) - gamma0_ID, & !< accumulated shear at last switch of stress sense (at current switch?) - accshear_ID, & - shearrate_ID, & - resolvedstress_ID - end enum - type :: tParameters real(pReal) :: & gdot0, & !< reference shear strain rate for slip n, & !< stress exponent for slip aTolResistance, & aTolShear - real(pReal), allocatable, dimension(:) :: & + real(pReal), allocatable, dimension(:) :: & crss0, & !< initial critical shear stress for slip theta0, & !< initial hardening rate of forward stress for each slip theta1, & !< asymptotic hardening rate of forward stress for each slip @@ -35,21 +22,21 @@ submodule(constitutive) plastic_kinehardening tau1, & tau1_b, & nonSchmidCoeff - real(pReal), allocatable, dimension(:,:) :: & + real(pReal), allocatable, dimension(:,:) :: & interaction_slipslip !< slip resistance from slip activity - real(pReal), allocatable, dimension(:,:,:) :: & + real(pReal), allocatable, dimension(:,:,:) :: & Schmid, & nonSchmid_pos, & nonSchmid_neg integer :: & totalNslip, & !< total number of active slip system of_debug = 0 - integer, allocatable, dimension(:) :: & + integer, allocatable, dimension(:) :: & Nslip !< number of active slip systems for each family - integer(kind(undefined_ID)), allocatable, dimension(:) :: & - outputID !< ID of each post result output + character(len=pStringLen), allocatable, dimension(:) :: & + output end type tParameters - + type :: tKinehardeningState real(pReal), pointer, dimension(:,:) :: & !< vectors along NipcMyInstance crss, & !< critical resolved stress @@ -59,7 +46,7 @@ submodule(constitutive) plastic_kinehardening gamma0, & !< accumulated shear at last switch of stress sense accshear !< accumulated (absolute) shear end type tKinehardeningState - + !-------------------------------------------------------------------------------------------------- ! containers for parameters and state type(tParameters), allocatable, dimension(:) :: param @@ -72,37 +59,32 @@ contains !-------------------------------------------------------------------------------------------------- -!> @brief module initialization +!> @brief Perform module initialization. !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- module subroutine plastic_kinehardening_init integer :: & Ninstance, & - p, i, o, & + p, o, & NipcMyPhase, & sizeState, sizeDeltaState, sizeDotState, & startIndex, endIndex - - integer(kind(undefined_ID)) :: & - outputID - + character(len=pStringLen) :: & extmsg = '' - character(len=pStringLen), dimension(:), allocatable :: & - outputs - - write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_KINEHARDENING_label//' init -+>>>' - + + write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_KINEHARDENING_label//' init -+>>>'; flush(6) + Ninstance = count(phase_plasticity == PLASTICITY_KINEHARDENING_ID) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - + allocate(param(Ninstance)) allocate(state(Ninstance)) allocate(dotState(Ninstance)) allocate(deltaState(Ninstance)) - + do p = 1, size(phase_plasticityInstance) if (phase_plasticity(p) /= PLASTICITY_KINEHARDENING_ID) cycle associate(prm => param(phase_plasticityInstance(p)), & @@ -110,22 +92,22 @@ module subroutine plastic_kinehardening_init dlt => deltaState(phase_plasticityInstance(p)), & stt => state(phase_plasticityInstance(p)),& config => config_phase(p)) - + #ifdef DEBUG if (p==material_phaseAt(debug_g,debug_e)) then - prm%of_debug = material_phasememberAt(debug_g,debug_i,debug_e) + prm%of_debug = material_phasememberAt(debug_g,debug_i,debug_e) endif #endif - + !-------------------------------------------------------------------------------------------------- ! optional parameters that need to be defined prm%aTolResistance = config%getFloat('atol_resistance',defaultVal=1.0_pReal) prm%aTolShear = config%getFloat('atol_shear', defaultVal=1.0e-6_pReal) - + ! sanity checks if (prm%aTolResistance <= 0.0_pReal) extmsg = trim(extmsg)//' aTolresistance' if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//' aTolShear' - + !-------------------------------------------------------------------------------------------------- ! slip related parameters prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) @@ -133,7 +115,7 @@ module subroutine plastic_kinehardening_init 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) @@ -146,7 +128,7 @@ module subroutine plastic_kinehardening_init prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(prm%Nslip, & config%getFloats('interaction_slipslip'), & config%getString('lattice_structure')) - + prm%crss0 = config%getFloats('crss0', requiredSize=size(prm%Nslip)) prm%tau1 = config%getFloats('tau1', requiredSize=size(prm%Nslip)) prm%tau1_b = config%getFloats('tau1_b', requiredSize=size(prm%Nslip)) @@ -154,10 +136,10 @@ module subroutine plastic_kinehardening_init prm%theta1 = config%getFloats('theta1', requiredSize=size(prm%Nslip)) prm%theta0_b = config%getFloats('theta0_b', requiredSize=size(prm%Nslip)) prm%theta1_b = config%getFloats('theta1_b', requiredSize=size(prm%Nslip)) - - prm%gdot0 = config%getFloat('gdot0') - prm%n = config%getFloat('n_slip') - + + prm%gdot0 = config%getFloat('gdot0') + prm%n = config%getFloat('n_slip') + ! expand: family => system prm%crss0 = math_expand(prm%crss0, prm%Nslip) prm%tau1 = math_expand(prm%tau1, prm%Nslip) @@ -166,9 +148,9 @@ module subroutine plastic_kinehardening_init prm%theta1 = math_expand(prm%theta1, prm%Nslip) prm%theta0_b = math_expand(prm%theta0_b,prm%Nslip) prm%theta1_b = math_expand(prm%theta1_b,prm%Nslip) - - - + + + !-------------------------------------------------------------------------------------------------- ! sanity checks if ( prm%gdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0' @@ -176,58 +158,29 @@ module subroutine plastic_kinehardening_init if (any(prm%crss0 <= 0.0_pReal)) extmsg = trim(extmsg)//' crss0' if (any(prm%tau1 <= 0.0_pReal)) extmsg = trim(extmsg)//' tau1' if (any(prm%tau1_b <= 0.0_pReal)) extmsg = trim(extmsg)//' tau1_b' - + !ToDo: Any sensible checks for theta? - + endif slipActive - + !-------------------------------------------------------------------------------------------------- ! exit if any parameter is out of range if (extmsg /= '') & call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_KINEHARDENING_label//')') - + !-------------------------------------------------------------------------------------------------- ! output pararameters - outputs = config%getStrings('(output)',defaultVal=emptyStringArray) - allocate(prm%outputID(0)) - do i=1, size(outputs) - outputID = undefined_ID - select case(outputs(i)) - - case ('resistance') - outputID = merge(crss_ID,undefined_ID,prm%totalNslip>0) - case ('accumulatedshear') - outputID = merge(accshear_ID,undefined_ID,prm%totalNslip>0) - case ('shearrate') - outputID = merge(shearrate_ID,undefined_ID,prm%totalNslip>0) - case ('resolvedstress') - outputID = merge(resolvedstress_ID,undefined_ID,prm%totalNslip>0) - case ('backstress') - outputID = merge(crss_back_ID,undefined_ID,prm%totalNslip>0) - case ('sense') - outputID = merge(sense_ID,undefined_ID,prm%totalNslip>0) - case ('chi0') - outputID = merge(chi0_ID,undefined_ID,prm%totalNslip>0) - case ('gamma0') - outputID = merge(gamma0_ID,undefined_ID,prm%totalNslip>0) - - end select - - if (outputID /= undefined_ID) then - prm%outputID = [prm%outputID , outputID] - endif - - enddo - + prm%output = config%getStrings('(output)',defaultVal=emptyStringArray) + !-------------------------------------------------------------------------------------------------- ! allocate state arrays NipcMyPhase = count(material_phaseAt == p) * discretization_nIP sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%totalNslip sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%totalNslip sizeState = sizeDotState + sizeDeltaState - + call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,sizeDeltaState) - + !-------------------------------------------------------------------------------------------------- ! locally defined state aliases and initialization of state0 and aTolState startIndex = 1 @@ -236,13 +189,13 @@ module subroutine plastic_kinehardening_init stt%crss = spread(prm%crss0, 2, NipcMyPhase) dot%crss => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance - + startIndex = endIndex + 1 endIndex = endIndex + prm%totalNslip stt%crss_back => plasticState(p)%state (startIndex:endIndex,:) dot%crss_back => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance - + startIndex = endIndex + 1 endIndex = endIndex + prm%totalNslip stt%accshear => plasticState(p)%state (startIndex:endIndex,:) @@ -250,34 +203,34 @@ module subroutine plastic_kinehardening_init plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear ! global alias plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) - + o = plasticState(p)%offsetDeltaState startIndex = endIndex + 1 endIndex = endIndex + prm%totalNslip stt%sense => plasticState(p)%state (startIndex :endIndex ,:) dlt%sense => plasticState(p)%deltaState(startIndex-o:endIndex-o,:) - + startIndex = endIndex + 1 endIndex = endIndex + prm%totalNslip stt%chi0 => plasticState(p)%state (startIndex :endIndex ,:) dlt%chi0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:) - + startIndex = endIndex + 1 endIndex = endIndex + prm%totalNslip stt%gamma0 => plasticState(p)%state (startIndex :endIndex ,:) dlt%gamma0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:) - + plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally - + end associate - + enddo end subroutine plastic_kinehardening_init !-------------------------------------------------------------------------------------------------- -!> @brief calculates plastic velocity gradient and its tangent +!> @brief Calculate plastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) @@ -285,26 +238,26 @@ pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,insta Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp !< derivative of Lp with respect to the Mandel stress - + real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & instance, & of - + integer :: & i,k,l,m,n real(pReal), dimension(param(instance)%totalNslip) :: & gdot_pos,gdot_neg, & dgdot_dtau_pos,dgdot_dtau_neg - + Lp = 0.0_pReal dLp_dMp = 0.0_pReal - + associate(prm => param(instance)) - + call kinetics(Mp,instance,of,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) - + do i = 1, prm%totalNslip 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) & @@ -312,14 +265,14 @@ pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,insta + 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) enddo - + end associate end subroutine plastic_kinehardening_LpAndItsTangent !-------------------------------------------------------------------------------------------------- -!> @brief calculates the rate of change of microstructure +!> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- module subroutine plastic_kinehardening_dotState(Mp,instance,of) @@ -328,40 +281,40 @@ module subroutine plastic_kinehardening_dotState(Mp,instance,of) integer, intent(in) :: & instance, & of - + real(pReal) :: & sumGamma real(pReal), dimension(param(instance)%totalNslip) :: & gdot_pos,gdot_neg - - + + associate(prm => param(instance), stt => state(instance), dot => dotState(instance)) - + call kinetics(Mp,instance,of,gdot_pos,gdot_neg) dot%accshear(:,of) = abs(gdot_pos+gdot_neg) sumGamma = sum(stt%accshear(:,of)) - - + + dot%crss(:,of) = matmul(prm%interaction_SlipSlip,dot%accshear(:,of)) & * ( prm%theta1 & + (prm%theta0 - prm%theta1 + prm%theta0*prm%theta1*sumGamma/prm%tau1) & * exp(-sumGamma*prm%theta0/prm%tau1) & ) - + dot%crss_back(:,of) = stt%sense(:,of)*dot%accshear(:,of) * & ( prm%theta1_b + & (prm%theta0_b - prm%theta1_b & + prm%theta0_b*prm%theta1_b/(prm%tau1_b+stt%chi0(:,of))*(stt%accshear(:,of)-stt%gamma0(:,of))& ) *exp(-(stt%accshear(:,of)-stt%gamma0(:,of)) *prm%theta0_b/(prm%tau1_b+stt%chi0(:,of))) & ) - + end associate end subroutine plastic_kinehardening_dotState !-------------------------------------------------------------------------------------------------- -!> @brief calculates (instantaneous) incremental change of microstructure +!> @brief Calculate (instantaneous) incremental change of microstructure. !-------------------------------------------------------------------------------------------------- module subroutine plastic_kinehardening_deltaState(Mp,instance,of) @@ -370,18 +323,18 @@ module subroutine plastic_kinehardening_deltaState(Mp,instance,of) integer, intent(in) :: & instance, & of - + real(pReal), dimension(param(instance)%totalNslip) :: & gdot_pos,gdot_neg, & sense - + associate(prm => param(instance), stt => state(instance), dlt => deltaState(instance)) - + call kinetics(Mp,instance,of,gdot_pos,gdot_neg) sense = merge(state(instance)%sense(:,of), & ! keep existing... sign(1.0_pReal,gdot_pos+gdot_neg), & ! ...or have a defined dEq0(gdot_pos+gdot_neg,1e-10_pReal)) ! current sense of shear direction - + #ifdef DEBUG if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0 & .and. (of == prm%of_debug & @@ -390,7 +343,7 @@ module subroutine plastic_kinehardening_deltaState(Mp,instance,of) write(6,*) sense,state(instance)%sense(:,of) endif #endif - + !-------------------------------------------------------------------------------------------------- ! switch in sense of shear? where(dNeq(sense,stt%sense(:,of),0.1_pReal)) @@ -402,45 +355,43 @@ module subroutine plastic_kinehardening_deltaState(Mp,instance,of) dlt%chi0 (:,of) = 0.0_pReal dlt%gamma0(:,of) = 0.0_pReal end where - + end associate end subroutine plastic_kinehardening_deltaState !-------------------------------------------------------------------------------------------------- -!> @brief writes results to HDF5 output file +!> @brief Write results to HDF5 output file. !-------------------------------------------------------------------------------------------------- module subroutine plastic_kinehardening_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group + integer :: o associate(prm => param(instance), stt => state(instance)) - outputsLoop: do o = 1,size(prm%outputID) - select case(prm%outputID(o)) - case (crss_ID) - call results_writeDataset(group,stt%crss,'xi_sl', & - 'resistance against plastic slip','Pa') - - case(crss_back_ID) - call results_writeDataset(group,stt%crss_back,'tau_back', & - 'back stress against plastic slip','Pa') - - case (sense_ID) - call results_writeDataset(group,stt%sense,'sense_of_shear','tbd','1') - - case (chi0_ID) - call results_writeDataset(group,stt%chi0,'chi0','tbd','Pa') - - case (gamma0_ID) - call results_writeDataset(group,stt%gamma0,'gamma0','tbd','1') - - case (accshear_ID) - call results_writeDataset(group,stt%accshear,'gamma_sl', & - 'plastic shear','1') - + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) + case('resistance') + if(prm%totalNslip>0) call results_writeDataset(group,stt%crss,'xi_sl', & + 'resistance against plastic slip','Pa') + case('backstress') ! ToDo: should be 'tau_back' + if(prm%totalNslip>0) call results_writeDataset(group,stt%crss_back,'tau_back', & + 'back stress against plastic slip','Pa') + case ('sense') + if(prm%totalNslip>0) call results_writeDataset(group,stt%sense,'sense_of_shear', & + 'tbd','1') + case ('chi0') + if(prm%totalNslip>0) call results_writeDataset(group,stt%chi0,'chi0', & + 'tbd','Pa') + case ('gamma0') + if(prm%totalNslip>0) call results_writeDataset(group,stt%gamma0,'gamma0', & + 'tbd','1') + case ('accumulatedshear') + if(prm%totalNslip>0) call results_writeDataset(group,stt%accshear,'gamma_sl', & + 'plastic shear','1') end select enddo outputsLoop end associate @@ -449,10 +400,11 @@ end subroutine plastic_kinehardening_results !-------------------------------------------------------------------------------------------------- -!> @brief calculates shear rates on slip systems and derivatives with respect to resolved stress -!> @details: Shear rates are calculated only optionally. +!> @brief Calculate shear rates on slip systems and their derivatives with respect to resolved +! stress. +!> @details: Derivatives are calculated only optionally. ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to -! have the optional arguments at the end +! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- pure subroutine kinetics(Mp,instance,of, & gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) @@ -462,44 +414,44 @@ pure subroutine kinetics(Mp,instance,of, & integer, intent(in) :: & instance, & of - + real(pReal), intent(out), dimension(param(instance)%totalNslip) :: & gdot_pos, & gdot_neg real(pReal), intent(out), optional, dimension(param(instance)%totalNslip) :: & dgdot_dtau_pos, & dgdot_dtau_neg - + real(pReal), dimension(param(instance)%totalNslip) :: & tau_pos, & tau_neg integer :: i logical :: nonSchmidActive - + associate(prm => param(instance), stt => state(instance)) - + nonSchmidActive = size(prm%nonSchmidCoeff) > 0 - + do i = 1, prm%totalNslip tau_pos(i) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i)) - stt%crss_back(i,of) tau_neg(i) = merge(math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i)) - stt%crss_back(i,of), & 0.0_pReal, nonSchmidActive) enddo - + where(dNeq0(tau_pos)) gdot_pos = prm%gdot0 * merge(0.5_pReal,1.0_pReal, nonSchmidActive) & ! 1/2 if non-Schmid active * sign(abs(tau_pos/stt%crss(:,of))**prm%n, tau_pos) else where gdot_pos = 0.0_pReal end where - + where(dNeq0(tau_neg)) gdot_neg = prm%gdot0 * 0.5_pReal & ! only used if non-Schmid active, always 1/2 * sign(abs(tau_neg/stt%crss(:,of))**prm%n, tau_neg) else where gdot_neg = 0.0_pReal end where - + if (present(dgdot_dtau_pos)) then where(dNeq0(gdot_pos)) dgdot_dtau_pos = gdot_pos*prm%n/tau_pos diff --git a/src/constitutive_plastic_none.f90 b/src/constitutive_plastic_none.f90 index 578fb4149..4840c9c09 100644 --- a/src/constitutive_plastic_none.f90 +++ b/src/constitutive_plastic_none.f90 @@ -18,19 +18,19 @@ module subroutine plastic_none_init Ninstance, & p, & NipcMyPhase - - write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_NONE_label//' init -+>>>' - + + write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_NONE_label//' init -+>>>'; flush(6) + Ninstance = count(phase_plasticity == PLASTICITY_NONE_ID) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - + do p = 1, size(phase_plasticity) if (phase_plasticity(p) /= PLASTICITY_NONE_ID) cycle - + NipcMyPhase = count(material_phaseAt == p) * discretization_nIP call material_allocatePlasticState(p,NipcMyPhase,0,0,0) - + enddo end subroutine plastic_none_init diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index 3739ab669..47db19385 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -48,32 +48,6 @@ submodule(constitutive) plastic_nonlocal real(pReal), dimension(:,:,:,:,:,:), allocatable :: & compatibility !< slip system compatibility between me and my neighbors - enum, bind(c) - enumerator :: & - undefined_ID, & - rho_sgl_mob_edg_pos_ID, & - rho_sgl_mob_edg_neg_ID, & - rho_sgl_mob_scr_pos_ID, & - rho_sgl_mob_scr_neg_ID, & - rho_sgl_imm_edg_pos_ID, & - rho_sgl_imm_edg_neg_ID, & - rho_sgl_imm_scr_pos_ID, & - rho_sgl_imm_scr_neg_ID, & - rho_dip_edg_ID, & - rho_dip_scr_ID, & - rho_forest_ID, & - resolvedstress_back_ID, & - tau_pass_ID, & - rho_dot_sgl_ID, & - rho_dot_sgl_mobile_ID, & - rho_dot_dip_ID, & - v_edg_pos_ID, & - v_edg_neg_ID, & - v_scr_pos_ID, & - v_scr_neg_ID, & - gamma_ID - end enum - type :: tParameters !< container type for internal constitutive parameters real(pReal) :: & atomicVolume, & !< atomic volume @@ -135,14 +109,12 @@ submodule(constitutive) plastic_nonlocal integer, dimension(:) ,allocatable :: & Nslip,& colinearSystem !< colinear system to the active slip system (only valid for fcc!) - + character(len=pStringLen), allocatable, dimension(:) :: & + output logical :: & shortRangeStressCorrection, & !< flag indicating the use of the short range stress correction by a excess density gradient term probabilisticMultiplication - integer(kind(undefined_ID)), dimension(:), allocatable :: & - outputID !< ID of each post result output - end type tParameters type :: tNonlocalMicrostructure @@ -198,22 +170,19 @@ module subroutine plastic_nonlocal_init integer :: & sizeState, sizeDotState,sizeDependentState, sizeDeltaState, & maxNinstances, & - p, i, & + p, & l, & s1, s2, & s, & t, & c - integer(kind(undefined_ID)) :: & - outputID character(len=pStringLen) :: & extmsg = '', & structure - character(len=pStringLen), dimension(:), allocatable :: outputs integer :: NofMyPhase - write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_NONLOCAL_label//' init -+>>>' + write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_NONLOCAL_label//' init -+>>>'; flush(6) write(6,'(/,a)') ' Reuber et al., Acta Materialia 71:333–348, 2014' write(6,'(a)') ' https://doi.org/10.1016/j.actamat.2014.03.012' @@ -407,60 +376,7 @@ module subroutine plastic_nonlocal_init endif slipActive - outputs = config%getStrings('(output)',defaultVal=emptyStringArray) - allocate(prm%outputID(0)) - do i=1, size(outputs) - outputID = undefined_ID - select case(trim(outputs(i))) - case ('rho_sgl_mob_edg_pos') - outputID = merge(rho_sgl_mob_edg_pos_ID,undefined_ID,prm%totalNslip>0) - case ('rho_sgl_mob_edg_neg') - outputID = merge(rho_sgl_mob_edg_neg_ID,undefined_ID,prm%totalNslip>0) - case ('rho_sgl_mob_scr_pos') - outputID = merge(rho_sgl_mob_scr_pos_ID,undefined_ID,prm%totalNslip>0) - case ('rho_sgl_mob_scr_neg') - outputID = merge(rho_sgl_mob_scr_neg_ID,undefined_ID,prm%totalNslip>0) - case ('rho_sgl_imm_edg_pos') - outputID = merge(rho_sgl_imm_edg_pos_ID,undefined_ID,prm%totalNslip>0) - case ('rho_sgl_imm_edg_neg') - outputID = merge(rho_sgl_imm_edg_neg_ID,undefined_ID,prm%totalNslip>0) - case ('rho_sgl_imm_scr_pos') - outputID = merge(rho_sgl_imm_scr_pos_ID,undefined_ID,prm%totalNslip>0) - case ('rho_sgl_imm_scr_neg') - outputID = merge(rho_sgl_imm_scr_neg_ID,undefined_ID,prm%totalNslip>0) - case ('rho_dip_edg') - outputID = merge(rho_dip_edg_ID,undefined_ID,prm%totalNslip>0) - case ('rho_dip_scr') - outputID = merge(rho_dip_scr_ID,undefined_ID,prm%totalNslip>0) - case ('rho_forest') - outputID = merge(rho_forest_ID,undefined_ID,prm%totalNslip>0) - case ('resolvedstress_back') - outputID = merge(resolvedstress_back_ID,undefined_ID,prm%totalNslip>0) - case ('tau_pass') - outputID = merge(tau_pass_ID,undefined_ID,prm%totalNslip>0) - case ('rho_dot_sgl') - outputID = merge(rho_dot_sgl_ID,undefined_ID,prm%totalNslip>0) - case ('rho_dot_sgl_mobile') - outputID = merge(rho_dot_sgl_mobile_ID,undefined_ID,prm%totalNslip>0) - case ('rho_dot_dip') - outputID = merge(rho_dot_dip_ID,undefined_ID,prm%totalNslip>0) - case ('v_edg_pos') - outputID = merge(v_edg_pos_ID,undefined_ID,prm%totalNslip>0) - case ('v_edg_neg') - outputID = merge(v_edg_neg_ID,undefined_ID,prm%totalNslip>0) - case ('v_scr_pos') - outputID = merge(v_scr_pos_ID,undefined_ID,prm%totalNslip>0) - case ('v_scr_neg') - outputID = merge(v_scr_neg_ID,undefined_ID,prm%totalNslip>0) - case ('gamma') - outputID = merge(gamma_ID,undefined_ID,prm%totalNslip>0) - end select - - if (outputID /= undefined_ID) then - prm%outputID = [prm%outputID , outputID] - endif - - enddo + prm%output = config%getStrings('(output)',defaultVal=emptyStringArray) !-------------------------------------------------------------------------------------------------- ! allocate state arrays @@ -561,7 +477,7 @@ module subroutine plastic_nonlocal_init stt%v_scr_neg => plasticState(p)%state (15*prm%totalNslip + 1:16*prm%totalNslip ,1:NofMyPhase) allocate(dst%tau_pass(prm%totalNslip,NofMyPhase),source=0.0_pReal) - allocate(dst%tau_Back(prm%totalNslip,NofMyPhase), source=0.0_pReal) + allocate(dst%tau_back(prm%totalNslip,NofMyPhase),source=0.0_pReal) end associate if (NofMyPhase > 0) call stateInit(p,NofMyPhase) @@ -1968,62 +1884,63 @@ module subroutine plastic_nonlocal_results(instance,group) integer, intent(in) :: instance character(len=*),intent(in) :: group + integer :: o associate(prm => param(instance),dst => microstructure(instance),stt=>state(instance)) - outputsLoop: do o = 1,size(prm%outputID) - select case(prm%outputID(o)) - case (rho_sgl_mob_edg_pos_ID) - call results_writeDataset(group,stt%rho_sgl_mob_edg_pos, 'rho_sgl_mob_edg_pos', & - 'positive mobile edge density','1/m²') - case (rho_sgl_imm_edg_pos_ID) - call results_writeDataset(group,stt%rho_sgl_imm_edg_pos, 'rho_sgl_imm_edg_pos',& - 'positive immobile edge density','1/m²') - case (rho_sgl_mob_edg_neg_ID) - call results_writeDataset(group,stt%rho_sgl_mob_edg_neg, 'rho_sgl_mob_edg_neg',& - 'negative mobile edge density','1/m²') - case (rho_sgl_imm_edg_neg_ID) - call results_writeDataset(group,stt%rho_sgl_imm_edg_neg, 'rho_sgl_imm_edg_neg',& - 'negative immobile edge density','1/m²') - case (rho_dip_edg_ID) - call results_writeDataset(group,stt%rho_dip_edg, 'rho_dip_edg',& - 'edge dipole density','1/m²') - case (rho_sgl_mob_scr_pos_ID) - call results_writeDataset(group,stt%rho_sgl_mob_scr_pos, 'rho_sgl_mob_scr_pos',& - 'positive mobile screw density','1/m²') - case (rho_sgl_imm_scr_pos_ID) - call results_writeDataset(group,stt%rho_sgl_imm_scr_pos, 'rho_sgl_imm_scr_pos',& - 'positive immobile screw density','1/m²') - case (rho_sgl_mob_scr_neg_ID) - call results_writeDataset(group,stt%rho_sgl_mob_scr_neg, 'rho_sgl_mob_scr_neg',& - 'negative mobile screw density','1/m²') - case (rho_sgl_imm_scr_neg_ID) - call results_writeDataset(group,stt%rho_sgl_imm_scr_neg, 'rho_sgl_imm_scr_neg',& - 'negative immobile screw density','1/m²') - case (rho_dip_scr_ID) - call results_writeDataset(group,stt%rho_dip_scr, 'rho_dip_scr',& - 'screw dipole density','1/m²') - case (rho_forest_ID) - call results_writeDataset(group,stt%rho_forest, 'rho_forest',& - 'forest density','1/m²') - case (v_edg_pos_ID) - call results_writeDataset(group,stt%v_edg_pos, 'v_edg_pos',& - 'positive edge velocity','m/s') - case (v_edg_neg_ID) - call results_writeDataset(group,stt%v_edg_neg, 'v_edg_neg',& - 'negative edge velocity','m/s') - case (v_scr_pos_ID) - call results_writeDataset(group,stt%v_scr_pos, 'v_scr_pos',& - 'positive srew velocity','m/s') - case (v_scr_neg_ID) - call results_writeDataset(group,stt%v_scr_neg, 'v_scr_neg',& - 'negative screw velocity','m/s') - case(gamma_ID) - call results_writeDataset(group,stt%gamma,'gamma',& - 'plastic shear','1') - case (tau_pass_ID) - call results_writeDataset(group,dst%tau_pass,'tau_pass',& - 'passing stress for slip','Pa') + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) + case('rho_sgl_mob_edg_pos') + if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_sgl_mob_edg_pos, 'rho_sgl_mob_edg_pos', & + 'positive mobile edge density','1/m²') + case('rho_sgl_imm_edg_pos') + if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_sgl_imm_edg_pos, 'rho_sgl_imm_edg_pos',& + 'positive immobile edge density','1/m²') + case('rho_sgl_mob_edg_neg') + if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_sgl_mob_edg_neg, 'rho_sgl_mob_edg_neg',& + 'negative mobile edge density','1/m²') + case('rho_sgl_imm_edg_neg') + if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_sgl_imm_edg_neg, 'rho_sgl_imm_edg_neg',& + 'negative immobile edge density','1/m²') + case('rho_dip_edg') + if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_dip_edg, 'rho_dip_edg',& + 'edge dipole density','1/m²') + case('rho_sgl_mob_scr_pos') + if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_sgl_mob_scr_pos, 'rho_sgl_mob_scr_pos',& + 'positive mobile screw density','1/m²') + case('rho_sgl_imm_scr_pos') + if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_sgl_imm_scr_pos, 'rho_sgl_imm_scr_pos',& + 'positive immobile screw density','1/m²') + case('rho_sgl_mob_scr_neg') + if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_sgl_mob_scr_neg, 'rho_sgl_mob_scr_neg',& + 'negative mobile screw density','1/m²') + case('rho_sgl_imm_scr_neg') + if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_sgl_imm_scr_neg, 'rho_sgl_imm_scr_neg',& + 'negative immobile screw density','1/m²') + case('rho_dip_scr') + if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_dip_scr, 'rho_dip_scr',& + 'screw dipole density','1/m²') + case('rho_forest') + if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_forest, 'rho_forest',& + 'forest density','1/m²') + case('v_edg_pos') + if(prm%totalNslip>0) call results_writeDataset(group,stt%v_edg_pos, 'v_edg_pos',& + 'positive edge velocity','m/s') + case('v_edg_neg') + if(prm%totalNslip>0) call results_writeDataset(group,stt%v_edg_neg, 'v_edg_neg',& + 'negative edge velocity','m/s') + case('v_scr_pos') + if(prm%totalNslip>0) call results_writeDataset(group,stt%v_scr_pos, 'v_scr_pos',& + 'positive srew velocity','m/s') + case('v_scr_neg') + if(prm%totalNslip>0) call results_writeDataset(group,stt%v_scr_neg, 'v_scr_neg',& + 'negative screw velocity','m/s') + case('gamma') + if(prm%totalNslip>0) call results_writeDataset(group,stt%gamma,'gamma',& + 'plastic shear','1') + case('tau_pass') + if(prm%totalNslip>0) call results_writeDataset(group,dst%tau_pass,'tau_pass',& + 'passing stress for slip','Pa') end select enddo outputsLoop end associate diff --git a/src/constitutive_plastic_phenopowerlaw.f90 b/src/constitutive_plastic_phenopowerlaw.f90 index 165ad0081..6dfa3fb16 100644 --- a/src/constitutive_plastic_phenopowerlaw.f90 +++ b/src/constitutive_plastic_phenopowerlaw.f90 @@ -6,19 +6,6 @@ !-------------------------------------------------------------------------------------------------- submodule(constitutive) plastic_phenopowerlaw - enum, bind(c) - enumerator :: & - undefined_ID, & - resistance_slip_ID, & - accumulatedshear_slip_ID, & - shearrate_slip_ID, & - resolvedstress_slip_ID, & - resistance_twin_ID, & - accumulatedshear_twin_ID, & - shearrate_twin_ID, & - resolvedstress_twin_ID - end enum - type :: tParameters real(pReal) :: & gdot0_slip, & !< reference shear strain rate for slip @@ -37,19 +24,19 @@ submodule(constitutive) plastic_phenopowerlaw aTolResistance, & !< absolute tolerance for integration of xi aTolShear, & !< absolute tolerance for integration of gamma aTolTwinfrac !< absolute tolerance for integration of f - real(pReal), allocatable, dimension(:) :: & + real(pReal), allocatable, dimension(:) :: & xi_slip_0, & !< initial critical shear stress for slip xi_twin_0, & !< initial critical shear stress for twin xi_slip_sat, & !< maximum critical shear stress for slip nonSchmidCoeff, & H_int, & !< per family hardening activity (optional) gamma_twin_char !< characteristic shear for twins - real(pReal), allocatable, dimension(:,:) :: & + real(pReal), allocatable, dimension(:,:) :: & interaction_SlipSlip, & !< slip resistance from slip activity interaction_SlipTwin, & !< slip resistance from twin activity interaction_TwinSlip, & !< twin resistance from slip activity interaction_TwinTwin !< twin resistance from twin activity - real(pReal), allocatable, dimension(:,:,:) :: & + real(pReal), allocatable, dimension(:,:,:) :: & Schmid_slip, & Schmid_twin, & nonSchmid_pos, & @@ -57,13 +44,13 @@ submodule(constitutive) plastic_phenopowerlaw integer :: & totalNslip, & !< total number of active slip system totalNtwin !< total number of active twin systems - integer, allocatable, dimension(:) :: & + integer, allocatable, dimension(:) :: & Nslip, & !< number of active slip systems for each family Ntwin !< number of active twin systems for each family - integer(kind(undefined_ID)), allocatable, dimension(:) :: & - outputID !< ID of each post result output + character(len=pStringLen), allocatable, dimension(:) :: & + output end type tParameters - + type :: tPhenopowerlawState real(pReal), pointer, dimension(:,:) :: & xi_slip, & @@ -71,7 +58,7 @@ submodule(constitutive) plastic_phenopowerlaw gamma_slip, & gamma_twin end type tPhenopowerlawState - + !-------------------------------------------------------------------------------------------------- ! containers for parameters and state type(tParameters), allocatable, dimension(:) :: param @@ -83,7 +70,7 @@ contains !-------------------------------------------------------------------------------------------------- -!> @brief module initialization +!> @brief Perform module initialization. !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- module subroutine plastic_phenopowerlaw_init @@ -91,51 +78,46 @@ module subroutine plastic_phenopowerlaw_init integer :: & Ninstance, & p, i, & - NipcMyPhase, outputSize, & + NipcMyPhase, & sizeState, sizeDotState, & startIndex, endIndex - - integer(kind(undefined_ID)) :: & - outputID - + character(len=pStringLen) :: & extmsg = '' - character(len=pStringLen), dimension(:), allocatable :: & - outputs - - write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_PHENOPOWERLAW_label//' init -+>>>' - + + write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_PHENOPOWERLAW_label//' init -+>>>'; flush(6) + Ninstance = count(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - + allocate(param(Ninstance)) allocate(state(Ninstance)) allocate(dotState(Ninstance)) - + do p = 1, size(phase_plasticity) if (phase_plasticity(p) /= PLASTICITY_PHENOPOWERLAW_ID) cycle associate(prm => param(phase_plasticityInstance(p)), & dot => dotState(phase_plasticityInstance(p)), & stt => state(phase_plasticityInstance(p)), & config => config_phase(p)) - + !-------------------------------------------------------------------------------------------------- ! optional parameters that need to be defined prm%c_1 = config%getFloat('twin_c',defaultVal=0.0_pReal) prm%c_2 = config%getFloat('twin_b',defaultVal=1.0_pReal) prm%c_3 = config%getFloat('twin_e',defaultVal=0.0_pReal) prm%c_4 = config%getFloat('twin_d',defaultVal=0.0_pReal) - + prm%aTolResistance = config%getFloat('atol_resistance',defaultVal=1.0_pReal) prm%aTolShear = config%getFloat('atol_shear', defaultVal=1.0e-6_pReal) prm%aTolTwinfrac = config%getFloat('atol_twinfrac', defaultVal=1.0e-6_pReal) - + ! sanity checks if (prm%aTolResistance <= 0.0_pReal) extmsg = trim(extmsg)//' aTolresistance' if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//' aTolShear' if (prm%aTolTwinfrac <= 0.0_pReal) extmsg = trim(extmsg)//' atoltwinfrac' - + !-------------------------------------------------------------------------------------------------- ! slip related parameters prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) @@ -143,7 +125,7 @@ module subroutine plastic_phenopowerlaw_init slipActive: if (prm%totalNslip > 0) then prm%Schmid_slip = 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) @@ -157,22 +139,22 @@ module subroutine plastic_phenopowerlaw_init prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(prm%Nslip, & config%getFloats('interaction_slipslip'), & config%getString('lattice_structure')) - + prm%xi_slip_0 = config%getFloats('tau0_slip', requiredSize=size(prm%Nslip)) prm%xi_slip_sat = config%getFloats('tausat_slip', requiredSize=size(prm%Nslip)) prm%H_int = config%getFloats('h_int', requiredSize=size(prm%Nslip), & defaultVal=[(0.0_pReal,i=1,size(prm%Nslip))]) - + prm%gdot0_slip = config%getFloat('gdot0_slip') prm%n_slip = config%getFloat('n_slip') prm%a_slip = config%getFloat('a_slip') prm%h0_SlipSlip = config%getFloat('h0_slipslip') - + ! expand: family => system prm%xi_slip_0 = math_expand(prm%xi_slip_0, prm%Nslip) prm%xi_slip_sat = math_expand(prm%xi_slip_sat,prm%Nslip) prm%H_int = math_expand(prm%H_int, prm%Nslip) - + ! sanity checks if ( prm%gdot0_slip <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0_slip' if ( prm%a_slip <= 0.0_pReal) extmsg = trim(extmsg)//' a_slip' @@ -183,7 +165,7 @@ module subroutine plastic_phenopowerlaw_init allocate(prm%interaction_SlipSlip(0,0)) allocate(prm%xi_slip_0(0)) endif slipActive - + !-------------------------------------------------------------------------------------------------- ! twin related parameters prm%Ntwin = config%getInts('ntwin', defaultVal=emptyIntArray) @@ -196,17 +178,17 @@ module subroutine plastic_phenopowerlaw_init config%getString('lattice_structure')) prm%gamma_twin_char = lattice_characteristicShear_twin(prm%Ntwin,config%getString('lattice_structure'),& config%getFloat('c/a')) - + prm%xi_twin_0 = config%getFloats('tau0_twin',requiredSize=size(prm%Ntwin)) - + prm%gdot0_twin = config%getFloat('gdot0_twin') prm%n_twin = config%getFloat('n_twin') prm%spr = config%getFloat('s_pr') prm%h0_TwinTwin = config%getFloat('h0_twintwin') - + ! expand: family => system prm%xi_twin_0 = math_expand(prm%xi_twin_0, prm%Ntwin) - + ! sanity checks if (prm%gdot0_twin <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0_twin' if (prm%n_twin <= 0.0_pReal) extmsg = trim(extmsg)//' n_twin' @@ -215,7 +197,7 @@ module subroutine plastic_phenopowerlaw_init allocate(prm%xi_twin_0(0)) allocate(prm%gamma_twin_char(0)) endif twinActive - + !-------------------------------------------------------------------------------------------------- ! slip-twin related parameters slipAndTwinActive: if (prm%totalNslip > 0 .and. prm%totalNtwin > 0) then @@ -231,63 +213,25 @@ module subroutine plastic_phenopowerlaw_init allocate(prm%interaction_TwinSlip(prm%TotalNtwin,prm%TotalNslip)) ! at least one dimension is 0 prm%h0_TwinSlip = 0.0_pReal endif slipAndTwinActive - + !-------------------------------------------------------------------------------------------------- ! exit if any parameter is out of range if (extmsg /= '') & call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_PHENOPOWERLAW_label//')') - + !-------------------------------------------------------------------------------------------------- ! output pararameters - outputs = config%getStrings('(output)',defaultVal=emptyStringArray) - allocate(prm%outputID(0)) - do i=1, size(outputs) - outputID = undefined_ID - select case(outputs(i)) - - case ('resistance_slip') - outputID = merge(resistance_slip_ID,undefined_ID,prm%totalNslip>0) - outputSize = prm%totalNslip - case ('accumulatedshear_slip') - outputID = merge(accumulatedshear_slip_ID,undefined_ID,prm%totalNslip>0) - outputSize = prm%totalNslip - case ('shearrate_slip') - outputID = merge(shearrate_slip_ID,undefined_ID,prm%totalNslip>0) - outputSize = prm%totalNslip - case ('resolvedstress_slip') - outputID = merge(resolvedstress_slip_ID,undefined_ID,prm%totalNslip>0) - outputSize = prm%totalNslip - - case ('resistance_twin') - outputID = merge(resistance_twin_ID,undefined_ID,prm%totalNtwin>0) - outputSize = prm%totalNtwin - case ('accumulatedshear_twin') - outputID = merge(accumulatedshear_twin_ID,undefined_ID,prm%totalNtwin>0) - outputSize = prm%totalNtwin - case ('shearrate_twin') - outputID = merge(shearrate_twin_ID,undefined_ID,prm%totalNtwin>0) - outputSize = prm%totalNtwin - case ('resolvedstress_twin') - outputID = merge(resolvedstress_twin_ID,undefined_ID,prm%totalNtwin>0) - outputSize = prm%totalNtwin - - end select - - if (outputID /= undefined_ID) then - prm%outputID = [prm%outputID, outputID] - endif - - enddo - + prm%output = config%getStrings('(output)',defaultVal=emptyStringArray) + !-------------------------------------------------------------------------------------------------- ! allocate state arrays NipcMyPhase = count(material_phaseAt == p) * discretization_nIP sizeDotState = size(['tau_slip ','gamma_slip']) * prm%totalNslip & + size(['tau_twin ','gamma_twin']) * prm%totalNtwin sizeState = sizeDotState - + call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0) - + !-------------------------------------------------------------------------------------------------- ! locally defined state aliases and initialization of state0 and aTolState startIndex = 1 @@ -296,14 +240,14 @@ module subroutine plastic_phenopowerlaw_init stt%xi_slip = spread(prm%xi_slip_0, 2, NipcMyPhase) dot%xi_slip => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance - + startIndex = endIndex + 1 endIndex = endIndex + prm%totalNtwin stt%xi_twin => plasticState(p)%state (startIndex:endIndex,:) stt%xi_twin = spread(prm%xi_twin_0, 2, NipcMyPhase) dot%xi_twin => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance - + startIndex = endIndex + 1 endIndex = endIndex + prm%totalNslip stt%gamma_slip => plasticState(p)%state (startIndex:endIndex,:) @@ -317,18 +261,18 @@ module subroutine plastic_phenopowerlaw_init stt%gamma_twin => plasticState(p)%state (startIndex:endIndex,:) dot%gamma_twin => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear - + plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally - + end associate - + enddo end subroutine plastic_phenopowerlaw_init !-------------------------------------------------------------------------------------------------- -!> @brief calculates plastic velocity gradient and its tangent +!> @brief Calculate plastic velocity gradient and its tangent. !> @details asummes that deformation by dislocation glide affects twinned and untwinned volume ! equally (Taylor assumption). Twinning happens only in untwinned volume !-------------------------------------------------------------------------------------------------- @@ -338,13 +282,13 @@ pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,insta Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp !< derivative of Lp with respect to the Mandel stress - + real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & instance, & of - + integer :: & i,k,l,m,n real(pReal), dimension(param(instance)%totalNslip) :: & @@ -352,12 +296,12 @@ pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,insta dgdot_dtauslip_pos,dgdot_dtauslip_neg real(pReal), dimension(param(instance)%totalNtwin) :: & gdot_twin,dgdot_dtautwin - + Lp = 0.0_pReal dLp_dMp = 0.0_pReal - + associate(prm => param(instance)) - + call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg) slipSystems: do i = 1, prm%totalNslip Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%Schmid_slip(1:3,1:3,i) @@ -366,7 +310,7 @@ pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,insta + dgdot_dtauslip_pos(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_pos(m,n,i) & + dgdot_dtauslip_neg(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_neg(m,n,i) enddo slipSystems - + call kinetics_twin(Mp,instance,of,gdot_twin,dgdot_dtautwin) twinSystems: do i = 1, prm%totalNtwin Lp = Lp + gdot_twin(i)*prm%Schmid_twin(1:3,1:3,i) @@ -374,14 +318,14 @@ pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,insta dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + dgdot_dtautwin(i)*prm%Schmid_twin(k,l,i)*prm%Schmid_twin(m,n,i) enddo twinSystems - + end associate end subroutine plastic_phenopowerlaw_LpAndItsTangent !-------------------------------------------------------------------------------------------------- -!> @brief calculates the rate of change of microstructure +!> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) @@ -390,7 +334,7 @@ module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) integer, intent(in) :: & instance, & of - + real(pReal) :: & c_SlipSlip,c_TwinSlip,c_TwinTwin, & xi_slip_sat_offset,& @@ -398,9 +342,9 @@ module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) real(pReal), dimension(param(instance)%totalNslip) :: & left_SlipSlip,right_SlipSlip, & gdot_slip_pos,gdot_slip_neg - + associate(prm => param(instance), stt => state(instance), dot => dotState(instance)) - + sumGamma = sum(stt%gamma_slip(:,of)) sumF = sum(stt%gamma_twin(:,of)/prm%gamma_twin_char) @@ -409,26 +353,26 @@ module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) c_SlipSlip = prm%h0_slipslip * (1.0_pReal + prm%c_1*sumF** prm%c_2) c_TwinSlip = prm%h0_TwinSlip * sumGamma**prm%c_3 c_TwinTwin = prm%h0_TwinTwin * sumF**prm%c_4 - + !-------------------------------------------------------------------------------------------------- ! calculate left and right vectors left_SlipSlip = 1.0_pReal + prm%H_int xi_slip_sat_offset = prm%spr*sqrt(sumF) right_SlipSlip = abs(1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset)) **prm%a_slip & * sign(1.0_pReal,1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset)) - + !-------------------------------------------------------------------------------------------------- ! shear rates call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg) dot%gamma_slip(:,of) = abs(gdot_slip_pos+gdot_slip_neg) call kinetics_twin(Mp,instance,of,dot%gamma_twin(:,of)) - + !-------------------------------------------------------------------------------------------------- ! hardening dot%xi_slip(:,of) = c_SlipSlip * left_SlipSlip * & matmul(prm%interaction_SlipSlip,dot%gamma_slip(:,of)*right_SlipSlip) & + matmul(prm%interaction_SlipTwin,dot%gamma_twin(:,of)) - + dot%xi_twin(:,of) = c_TwinSlip * matmul(prm%interaction_TwinSlip,dot%gamma_slip(:,of)) & + c_TwinTwin * matmul(prm%interaction_TwinTwin,dot%gamma_twin(:,of)) end associate @@ -437,33 +381,33 @@ end subroutine plastic_phenopowerlaw_dotState !-------------------------------------------------------------------------------------------------- -!> @brief writes results to HDF5 output file +!> @brief Write results to HDF5 output file. !-------------------------------------------------------------------------------------------------- module subroutine plastic_phenopowerlaw_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group - + integer :: o associate(prm => param(instance), stt => state(instance)) - outputsLoop: do o = 1,size(prm%outputID) - select case(prm%outputID(o)) + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) + + case('resistance_slip') + if(prm%totalNslip>0) call results_writeDataset(group,stt%xi_slip, 'xi_sl', & + 'resistance against plastic slip','Pa') + case('accumulatedshear_slip') + if(prm%totalNslip>0) call results_writeDataset(group,stt%gamma_slip,'gamma_sl', & + 'plastic shear','1') + + case('resistance_twin') + if(prm%totalNtwin>0) call results_writeDataset(group,stt%xi_twin, 'xi_tw', & + 'resistance against twinning','Pa') + case('accumulatedshear_twin') + if(prm%totalNtwin>0) call results_writeDataset(group,stt%gamma_twin,'gamma_tw', & + 'twinning shear','1') - case (resistance_slip_ID) - call results_writeDataset(group,stt%xi_slip, 'xi_sl', & - 'resistance against plastic slip','Pa') - case (accumulatedshear_slip_ID) - call results_writeDataset(group,stt%gamma_slip,'gamma_sl', & - 'plastic shear','1') - - case (resistance_twin_ID) - call results_writeDataset(group,stt%xi_twin, 'xi_tw', & - 'resistance against twinning','Pa') - case (accumulatedshear_twin_ID) - call results_writeDataset(group,stt%gamma_twin,'gamma_tw', & - 'twinning shear','1') - end select enddo outputsLoop end associate @@ -472,10 +416,11 @@ end subroutine plastic_phenopowerlaw_results !-------------------------------------------------------------------------------------------------- -!> @brief Shear rates on slip systems and their derivatives with respect to resolved stress +!> @brief Calculate shear rates on slip systems and their derivatives with respect to resolved +! stress. !> @details Derivatives are calculated only optionally. ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to -! have the optional arguments at the end +! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- pure subroutine kinetics_slip(Mp,instance,of, & gdot_slip_pos,gdot_slip_neg,dgdot_dtau_slip_pos,dgdot_dtau_slip_neg) @@ -485,44 +430,44 @@ pure subroutine kinetics_slip(Mp,instance,of, & integer, intent(in) :: & instance, & of - + real(pReal), intent(out), dimension(param(instance)%totalNslip) :: & gdot_slip_pos, & gdot_slip_neg real(pReal), intent(out), optional, dimension(param(instance)%totalNslip) :: & dgdot_dtau_slip_pos, & dgdot_dtau_slip_neg - + real(pReal), dimension(param(instance)%totalNslip) :: & tau_slip_pos, & tau_slip_neg integer :: i logical :: nonSchmidActive - + associate(prm => param(instance), stt => state(instance)) - + nonSchmidActive = size(prm%nonSchmidCoeff) > 0 - + do i = 1, prm%totalNslip tau_slip_pos(i) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i)) tau_slip_neg(i) = merge(math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i)), & 0.0_pReal, nonSchmidActive) enddo - + where(dNeq0(tau_slip_pos)) gdot_slip_pos = prm%gdot0_slip * merge(0.5_pReal,1.0_pReal, nonSchmidActive) & ! 1/2 if non-Schmid active * sign(abs(tau_slip_pos/stt%xi_slip(:,of))**prm%n_slip, tau_slip_pos) else where gdot_slip_pos = 0.0_pReal end where - + where(dNeq0(tau_slip_neg)) gdot_slip_neg = prm%gdot0_slip * 0.5_pReal & ! only used if non-Schmid active, always 1/2 * sign(abs(tau_slip_neg/stt%xi_slip(:,of))**prm%n_slip, tau_slip_neg) else where gdot_slip_neg = 0.0_pReal end where - + if (present(dgdot_dtau_slip_pos)) then where(dNeq0(gdot_slip_pos)) dgdot_dtau_slip_pos = gdot_slip_pos*prm%n_slip/tau_slip_pos @@ -543,9 +488,9 @@ end subroutine kinetics_slip !-------------------------------------------------------------------------------------------------- -!> @brief Shear rates on twin systems and their derivatives with respect to resolved stress. -! twinning is assumed to take place only in untwinned volume. -!> @details Derivates are calculated only optionally. +!> @brief Calculate shear rates on twin systems and their derivatives with respect to resolved +! stress. Twinning is assumed to take place only in untwinned volume. +!> @details Derivatives are calculated only optionally. ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- @@ -557,29 +502,29 @@ pure subroutine kinetics_twin(Mp,instance,of,& integer, intent(in) :: & instance, & of - + real(pReal), dimension(param(instance)%totalNtwin), intent(out) :: & gdot_twin real(pReal), dimension(param(instance)%totalNtwin), intent(out), optional :: & dgdot_dtau_twin - + real(pReal), dimension(param(instance)%totalNtwin) :: & tau_twin integer :: i - + associate(prm => param(instance), stt => state(instance)) - + do i = 1, prm%totalNtwin tau_twin(i) = math_mul33xx33(Mp,prm%Schmid_twin(1:3,1:3,i)) enddo - + where(tau_twin > 0.0_pReal) gdot_twin = (1.0_pReal-sum(stt%gamma_twin(:,of)/prm%gamma_twin_char)) & ! only twin in untwinned volume fraction * prm%gdot0_twin*(abs(tau_twin)/stt%xi_twin(:,of))**prm%n_twin else where gdot_twin = 0.0_pReal end where - + if (present(dgdot_dtau_twin)) then where(dNeq0(gdot_twin)) dgdot_dtau_twin = gdot_twin*prm%n_twin/tau_twin @@ -587,7 +532,7 @@ pure subroutine kinetics_twin(Mp,instance,of,& dgdot_dtau_twin = 0.0_pReal end where endif - + end associate end subroutine kinetics_twin diff --git a/src/homogenization_mech_RGC.f90 b/src/homogenization_mech_RGC.f90 index 9cb7efdc8..1f1d9d6ec 100644 --- a/src/homogenization_mech_RGC.f90 +++ b/src/homogenization_mech_RGC.f90 @@ -9,17 +9,6 @@ submodule(homogenization) homogenization_mech_RGC use rotations - enum, bind(c) - enumerator :: & - undefined_ID, & - constitutivework_ID, & - penaltyenergy_ID, & - volumediscrepancy_ID, & - averagerelaxrate_ID,& - maximumrelaxrate_ID,& - magnitudemismatch_ID - end enum - type :: tParameters integer, dimension(:), allocatable :: & Nconstituents @@ -31,8 +20,8 @@ submodule(homogenization) homogenization_mech_RGC angles integer :: & of_debug = 0 - integer(kind(undefined_ID)), dimension(:), allocatable :: & - outputID + character(len=pStringLen), allocatable, dimension(:) :: & + output end type tParameters type :: tRGCstate @@ -71,23 +60,18 @@ module subroutine mech_RGC_init integer :: & Ninstance, & - h, i, & + h, & NofMyHomog, & sizeState, nIntFaceTot - integer(kind(undefined_ID)) :: & - outputID - - character(len=pStringLen), dimension(:), allocatable :: & - outputs - write(6,'(/,a)') ' <<<+- homogenization_'//HOMOGENIZATION_RGC_label//' init -+>>>' + write(6,'(/,a)') ' <<<+- homogenization_'//HOMOGENIZATION_RGC_label//' init -+>>>'; flush(6) - write(6,'(/,a)') ' Tjahjanto et al., International Journal of Material Forming 2(1):939–942, 2009' - write(6,'(a)') ' https://doi.org/10.1007/s12289-009-0619-1' + write(6,'(/,a)') ' Tjahjanto et al., International Journal of Material Forming 2(1):939–942, 2009' + write(6,'(a)') ' https://doi.org/10.1007/s12289-009-0619-1' - write(6,'(/,a)') ' Tjahjanto et al., Modelling and Simulation in Materials Science and Engineering 18:015006, 2010' - write(6,'(a)') ' https://doi.org/10.1088/0965-0393/18/1/015006' + write(6,'(/,a)') ' Tjahjanto et al., Modelling and Simulation in Materials Science and Engineering 18:015006, 2010' + write(6,'(a)') ' https://doi.org/10.1088/0965-0393/18/1/015006' Ninstance = count(homogenization_type == HOMOGENIZATION_RGC_ID) if (iand(debug_level(debug_HOMOGENIZATION),debug_levelBasic) /= 0) & @@ -123,34 +107,8 @@ module subroutine mech_RGC_init prm%dAlpha = config%getFloats('grainsize', requiredSize=3) prm%angles = config%getFloats('clusterorientation',requiredSize=3) - outputs = config%getStrings('(output)',defaultVal=emptyStringArray) - allocate(prm%outputID(0)) - - do i=1, size(outputs) - outputID = undefined_ID - select case(outputs(i)) - - case('constitutivework') - outputID = constitutivework_ID - case('penaltyenergy') - outputID = penaltyenergy_ID - case('volumediscrepancy') - outputID = volumediscrepancy_ID - case('averagerelaxrate') - outputID = averagerelaxrate_ID - case('maximumrelaxrate') - outputID = maximumrelaxrate_ID - case('magnitudemismatch') - outputID = magnitudemismatch_ID - - end select - - if (outputID /= undefined_ID) then - prm%outputID = [prm%outputID , outputID] - endif - - enddo - + prm%output = config%getStrings('(output)',defaultVal=emptyStringArray) + NofMyHomog = count(material_homogenizationAt == h) nIntFaceTot = 3*( (prm%Nconstituents(1)-1)*prm%Nconstituents(2)*prm%Nconstituents(3) & + prm%Nconstituents(1)*(prm%Nconstituents(2)-1)*prm%Nconstituents(3) & @@ -934,26 +892,24 @@ module subroutine mech_RGC_results(instance,group) integer :: o associate(stt => state(instance), dst => dependentState(instance), prm => param(instance)) - - outputsLoop: do o = 1,size(prm%outputID) - select case(prm%outputID(o)) - - case (constitutivework_ID) + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) + case('constitutivework') call results_writeDataset(group,stt%work,'W',& 'work density','J/m³') - case (magnitudemismatch_ID) + case('magnitudemismatch') call results_writeDataset(group,dst%mismatch,'N',& 'average mismatch tensor','1') - case (penaltyenergy_ID) + case('penaltyenergy') call results_writeDataset(group,stt%penaltyEnergy,'R',& 'mismatch penalty density','J/m³') - case (volumediscrepancy_ID) + case('volumediscrepancy') call results_writeDataset(group,dst%volumeDiscrepancy,'Delta_V',& 'volume discrepancy','m³') - case (maximumrelaxrate_ID) + case('maximumrelaxrate') call results_writeDataset(group,dst%relaxationrate_max,'max_alpha_dot',& 'maximum relaxation rate','m/s') - case (averagerelaxrate_ID) + case('averagerelaxrate') call results_writeDataset(group,dst%relaxationrate_avg,'avg_alpha_dot',& 'average relaxation rate','m/s') end select diff --git a/src/thermal_adiabatic.f90 b/src/thermal_adiabatic.f90 index 5fa6b3281..2bef1c3ee 100644 --- a/src/thermal_adiabatic.f90 +++ b/src/thermal_adiabatic.f90 @@ -16,14 +16,9 @@ module thermal_adiabatic implicit none private - enum, bind(c) - enumerator :: undefined_ID, & - temperature_ID - end enum - type :: tParameters - integer(kind(undefined_ID)), dimension(:), allocatable :: & - outputID + character(len=pStringLen), allocatable, dimension(:) :: & + output end type tParameters type(tparameters), dimension(:), allocatable :: & @@ -46,10 +41,9 @@ contains !-------------------------------------------------------------------------------------------------- subroutine thermal_adiabatic_init - integer :: maxNinstance,o,h,NofMyHomog - character(len=pStringLen), dimension(:), allocatable :: outputs + integer :: maxNinstance,h,NofMyHomog - write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_ADIABATIC_label//' init -+>>>'; flush(6) + write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_ADIABATIC_label//' init -+>>>'; flush(6) maxNinstance = count(thermal_type == THERMAL_adiabatic_ID) if (maxNinstance == 0) return @@ -60,15 +54,7 @@ subroutine thermal_adiabatic_init if (thermal_type(h) /= THERMAL_adiabatic_ID) cycle associate(prm => param(thermal_typeInstance(h)),config => config_homogenization(h)) - outputs = config%getStrings('(output)',defaultVal=emptyStringArray) - allocate(prm%outputID(0)) - - do o=1, size(outputs) - select case(outputs(o)) - case('temperature') - prm%outputID = [prm%outputID, temperature_ID] - end select - enddo + prm%output = config%getStrings('(output)',defaultVal=emptyStringArray) NofMyHomog=count(material_homogenizationAt==h) thermalState(h)%sizeState = 1 @@ -76,7 +62,6 @@ subroutine thermal_adiabatic_init allocate(thermalState(h)%subState0(1,NofMyHomog), source=thermal_initialT(h)) allocate(thermalState(h)%state (1,NofMyHomog), source=thermal_initialT(h)) - nullify(thermalMapping(h)%p) thermalMapping(h)%p => material_homogenizationMemberAt deallocate(temperature(h)%p) temperature(h)%p => thermalState(h)%state(1,:) @@ -246,14 +231,13 @@ subroutine thermal_adiabatic_results(homog,group) integer, intent(in) :: homog character(len=*), intent(in) :: group + integer :: o associate(prm => param(damage_typeInstance(homog))) - - outputsLoop: do o = 1,size(prm%outputID) - select case(prm%outputID(o)) - - case (temperature_ID) + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) + case('temperature') ! ToDo: should be 'T' call results_writeDataset(group,temperature(homog)%p,'T',& 'temperature','K') end select diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90 index 9ec8082ec..dad3fecd8 100644 --- a/src/thermal_conduction.f90 +++ b/src/thermal_conduction.f90 @@ -15,15 +15,9 @@ module thermal_conduction implicit none private - enum, bind(c) - enumerator :: & - undefined_ID, & - temperature_ID - end enum - type :: tParameters - integer(kind(undefined_ID)), dimension(:), allocatable :: & - outputID + character(len=pStringLen), allocatable, dimension(:) :: & + output end type tParameters type(tparameters), dimension(:), allocatable :: & @@ -48,10 +42,9 @@ contains subroutine thermal_conduction_init - integer :: maxNinstance,o,NofMyHomog,h - character(len=pStringLen), dimension(:), allocatable :: outputs + integer :: maxNinstance,NofMyHomog,h - write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_CONDUCTION_label//' init -+>>>'; flush(6) + write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_CONDUCTION_label//' init -+>>>'; flush(6) maxNinstance = count(thermal_type == THERMAL_conduction_ID) if (maxNinstance == 0) return @@ -62,15 +55,7 @@ subroutine thermal_conduction_init if (thermal_type(h) /= THERMAL_conduction_ID) cycle associate(prm => param(thermal_typeInstance(h)),config => config_homogenization(h)) - outputs = config%getStrings('(output)',defaultVal=emptyStringArray) - allocate(prm%outputID(0)) - - do o=1, size(outputs) - select case(outputs(o)) - case('temperature') - prm%outputID = [prm%outputID, temperature_ID] - end select - enddo + prm%output = config%getStrings('(output)',defaultVal=emptyStringArray) NofMyHomog=count(material_homogenizationAt==h) thermalState(h)%sizeState = 0 @@ -78,7 +63,6 @@ subroutine thermal_conduction_init allocate(thermalState(h)%subState0(0,NofMyHomog)) allocate(thermalState(h)%state (0,NofMyHomog)) - nullify(thermalMapping(h)%p) thermalMapping(h)%p => material_homogenizationMemberAt deallocate(temperature (h)%p) allocate (temperature (h)%p(NofMyHomog), source=thermal_initialT(h)) @@ -259,14 +243,13 @@ subroutine thermal_conduction_results(homog,group) integer, intent(in) :: homog character(len=*), intent(in) :: group + integer :: o associate(prm => param(damage_typeInstance(homog))) - - outputsLoop: do o = 1,size(prm%outputID) - select case(prm%outputID(o)) - - case (temperature_ID) + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) + case('temperature') ! ToDo: should be 'T' call results_writeDataset(group,temperature(homog)%p,'T',& 'temperature','K') end select