diff --git a/src/kinematics_slipplane_opening.f90 b/src/kinematics_slipplane_opening.f90 index 86be20c9d..0fa80ece8 100644 --- a/src/kinematics_slipplane_opening.f90 +++ b/src/kinematics_slipplane_opening.f90 @@ -22,26 +22,14 @@ module kinematics_slipplane_opening sdot0, & n real(pReal), dimension(:), allocatable :: & - critDisp, & - critPlasticStrain - end type + critLoad + real(pReal), dimension(:,:), allocatable :: & + slip_direction, & + slip_normal, & + slip_transverse + end type tParameters -! Begin Deprecated - integer(pInt), dimension(:), allocatable, private :: & - kinematics_slipplane_opening_totalNslip !< total number of slip systems - - integer(pInt), dimension(:,:), allocatable, private :: & - kinematics_slipplane_opening_Nslip !< number of slip systems per family - - real(pReal), dimension(:), allocatable, private :: & - kinematics_slipplane_opening_sdot_0, & - kinematics_slipplane_opening_N - - real(pReal), dimension(:,:), allocatable, private :: & - kinematics_slipplane_opening_critPlasticStrain, & - kinematics_slipplane_opening_critLoad -! End Deprecated - + type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) public :: & kinematics_slipplane_opening_init, & kinematics_slipplane_opening_LiAndItsTangent @@ -54,11 +42,6 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine kinematics_slipplane_opening_init() -#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 - use, intrinsic :: iso_fortran_env, only: & - compiler_version, & - compiler_options -#endif use debug, only: & debug_level,& debug_constitutive,& @@ -66,29 +49,23 @@ subroutine kinematics_slipplane_opening_init() use config, only: & config_phase use IO, only: & - IO_warning, & - IO_error, & - IO_timeStamp + IO_error + use math, only: & + math_expand use material, only: & phase_kinematics, & KINEMATICS_slipplane_opening_label, & KINEMATICS_slipplane_opening_ID - use lattice, only: & - lattice_maxNslipFamily, & - lattice_NslipSystem + use lattice implicit none - integer(pInt), allocatable, dimension(:) :: tempInt - real(pReal), allocatable, dimension(:) :: tempFloat integer(pInt) :: maxNinstance,p,instance,kinematics write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_slipplane_opening_LABEL//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() -#include "compilation_info.f90" - maxNinstance = int(count(phase_kinematics == KINEMATICS_slipplane_opening_ID),pInt) - if (maxNinstance == 0_pInt) return + maxNinstance = count(phase_kinematics == KINEMATICS_slipplane_opening_ID) + if (maxNinstance == 0) return if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance @@ -97,38 +74,38 @@ subroutine kinematics_slipplane_opening_init() do p = 1_pInt, size(config_phase) kinematics_slipplane_opening_instance(p) = count(phase_kinematics(:,1:p) == kinematics_slipplane_opening_ID) ! ToDo: count correct? enddo - - allocate(kinematics_slipplane_opening_critLoad(lattice_maxNslipFamily,maxNinstance), source=0.0_pReal) - allocate(kinematics_slipplane_opening_critPlasticStrain(lattice_maxNslipFamily,maxNinstance),source=0.0_pReal) - allocate(kinematics_slipplane_opening_Nslip(lattice_maxNslipFamily,maxNinstance), source=0_pInt) - allocate(kinematics_slipplane_opening_totalNslip(maxNinstance), source=0_pInt) - allocate(kinematics_slipplane_opening_N(maxNinstance), source=0.0_pReal) - allocate(kinematics_slipplane_opening_sdot_0(maxNinstance), source=0.0_pReal) - + + allocate(param(maxNinstance)) + do p = 1_pInt, size(config_phase) if (all(phase_kinematics(:,p) /= KINEMATICS_slipplane_opening_ID)) cycle + associate(prm => param(kinematics_slipplane_opening_instance(p)), & + config => config_phase(p)) instance = kinematics_slipplane_opening_instance(p) - kinematics_slipplane_opening_sdot_0(instance) = config_phase(p)%getFloat('anisoductile_sdot0') - kinematics_slipplane_opening_N(instance) = config_phase(p)%getFloat('anisoductile_ratesensitivity') - tempInt = config_phase(p)%getInts('ncleavage') - kinematics_slipplane_opening_Nslip(1:size(tempInt),instance) = tempInt + prm%sdot0 = config_phase(p)%getFloat('anisoductile_sdot0') + prm%n = config_phase(p)%getFloat('anisoductile_ratesensitivity') + + prm%Nslip = config%getInts('nslip') - tempFloat = config_phase(p)%getFloats('anisoductile_criticalplasticstrain',requiredSize=size(tempInt)) - kinematics_slipplane_opening_critPlasticStrain(1:size(tempInt),instance) = tempFloat + prm%critLoad = config_phase(p)%getFloats('anisoductile_criticalload',requiredSize=size(prm%Nslip )) - tempFloat = config_phase(p)%getFloats('anisoductile_criticalload',requiredSize=size(tempInt)) - kinematics_slipplane_opening_critLoad(1:size(tempInt),instance) = tempFloat + prm%critLoad = math_expand(prm%critLoad, prm%Nslip) + +prm%slip_direction = lattice_slip_direction (prm%Nslip,config%getString('lattice_structure'),& + config%getFloat('c/a',defaultVal=0.0_pReal)) + prm%slip_normal = lattice_slip_normal (prm%Nslip,config%getString('lattice_structure'),& + config%getFloat('c/a',defaultVal=0.0_pReal)) + prm%slip_transverse = lattice_slip_transverse(prm%Nslip,config%getString('lattice_structure'),& + config%getFloat('c/a',defaultVal=0.0_pReal)) - kinematics_slipplane_opening_Nslip(1:lattice_maxNslipFamily,instance) = & - min(lattice_NslipSystem(1:lattice_maxNslipFamily,p),& ! limit active cleavage systems per family to min of available and requested - kinematics_slipplane_opening_Nslip(1:lattice_maxNslipFamily,instance)) - kinematics_slipplane_opening_totalNslip(instance) = sum(kinematics_slipplane_opening_Nslip(:,instance)) - if (kinematics_slipplane_opening_sdot_0(instance) <= 0.0_pReal) & - call IO_error(211_pInt,el=instance,ext_msg='sdot_0 ('//KINEMATICS_slipplane_opening_LABEL//')') - if (any(kinematics_slipplane_opening_critPlasticStrain(:,instance) < 0.0_pReal)) & - call IO_error(211_pInt,el=instance,ext_msg='criticaPlasticStrain ('//KINEMATICS_slipplane_opening_LABEL//')') - if (kinematics_slipplane_opening_N(instance) <= 0.0_pReal) & - call IO_error(211_pInt,el=instance,ext_msg='rate_sensitivity ('//KINEMATICS_slipplane_opening_LABEL//')') + ! if (kinematics_slipplane_opening_sdot_0(instance) <= 0.0_pReal) & + ! call IO_error(211_pInt,el=instance,ext_msg='sdot_0 ('//KINEMATICS_slipplane_opening_LABEL//')') + ! if (any(kinematics_slipplane_opening_critPlasticStrain(:,instance) < 0.0_pReal)) & + ! call IO_error(211_pInt,el=instance,ext_msg='criticaPlasticStrain ('//KINEMATICS_slipplane_opening_LABEL//')') + ! if (kinematics_slipplane_opening_N(instance) <= 0.0_pReal) & + ! call IO_error(211_pInt,el=instance,ext_msg='rate_sensitivity ('//KINEMATICS_slipplane_opening_LABEL//')') + + end associate enddo end subroutine kinematics_slipplane_opening_init @@ -140,23 +117,16 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, use prec, only: & tol_math_check use math, only: & - math_mul33xx33 - use lattice, only: & - lattice_maxNslipFamily, & - lattice_NslipSystem, & - lattice_sd, & - lattice_st, & - lattice_sn + math_mul33xx33, & + math_outer use material, only: & material_phase, & material_homog, & damage, & damageMapping - use math, only: & - math_tensorproduct33 - + implicit none - integer(pInt), intent(in) :: & + integer, intent(in) :: & ipc, & !< grain number ip, & !< integration point number el !< element number @@ -168,10 +138,10 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) real(pReal), dimension(3,3) :: & projection_d, projection_t, projection_n !< projection modes 3x3 tensor - integer(pInt) :: & + integer :: & instance, phase, & homog, damageOffset, & - f, i, index_myFamily, k, l, m, n + i, k, l, m, n real(pReal) :: & traction_d, traction_t, traction_n, traction_crit, & udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt @@ -181,66 +151,60 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, homog = material_homog(ip,el) damageOffset = damageMapping(homog)%p(ip,el) + associate(prm => param(instance)) Ld = 0.0_pReal dLd_dTstar = 0.0_pReal - do f = 1_pInt,lattice_maxNslipFamily - index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,phase)) ! at which index starts my family - do i = 1_pInt,kinematics_slipplane_opening_Nslip(f,instance) ! process each (active) slip system in family - projection_d = math_tensorproduct33(lattice_sd(1:3,index_myFamily+i,phase),& - lattice_sn(1:3,index_myFamily+i,phase)) - projection_t = math_tensorproduct33(lattice_st(1:3,index_myFamily+i,phase),& - lattice_sn(1:3,index_myFamily+i,phase)) - projection_n = math_tensorproduct33(lattice_sn(1:3,index_myFamily+i,phase),& - lattice_sn(1:3,index_myFamily+i,phase)) + do i = 1, prm%totalNslip + projection_d = math_outer(prm%slip_direction(1:3,i),prm%slip_normal(1:3,i)) + projection_t = math_outer(prm%slip_transverse(1:3,i),prm%slip_normal(1:3,i)) + projection_n = math_outer(prm%slip_normal(1:3,i),prm%slip_normal(1:3,i)) traction_d = math_mul33xx33(S,projection_d) traction_t = math_mul33xx33(S,projection_t) traction_n = math_mul33xx33(S,projection_n) - traction_crit = kinematics_slipplane_opening_critLoad(f,instance)* & - damage(homog)%p(damageOffset) ! degrading critical load carrying capacity by damage + traction_crit = prm%critLoad(i)* damage(homog)%p(damageOffset) ! degrading critical load carrying capacity by damage - udotd = & - sign(1.0_pReal,traction_d)* & - kinematics_slipplane_opening_sdot_0(instance)* & + udotd = sign(1.0_pReal,traction_d)* & + prm%sdot0* & (abs(traction_d)/traction_crit - & - abs(traction_d)/kinematics_slipplane_opening_critLoad(f,instance))**kinematics_slipplane_opening_N(instance) + abs(traction_d)/prm%critLoad(i))**prm%n if (abs(udotd) > tol_math_check) then Ld = Ld + udotd*projection_d - dudotd_dt = udotd*kinematics_slipplane_opening_N(instance)/traction_d + dudotd_dt = udotd*prm%n/traction_d forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + & dudotd_dt*projection_d(k,l)*projection_d(m,n) endif - udott = & - sign(1.0_pReal,traction_t)* & - kinematics_slipplane_opening_sdot_0(instance)* & + udott = sign(1.0_pReal,traction_t)* & + prm%sdot0* & (abs(traction_t)/traction_crit - & - abs(traction_t)/kinematics_slipplane_opening_critLoad(f,instance))**kinematics_slipplane_opening_N(instance) + abs(traction_t)/prm%critLoad(i))**prm%n if (abs(udott) > tol_math_check) then Ld = Ld + udott*projection_t - dudott_dt = udott*kinematics_slipplane_opening_N(instance)/traction_t + dudott_dt = udott*prm%n/traction_t forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + & dudott_dt*projection_t(k,l)*projection_t(m,n) endif udotn = & - kinematics_slipplane_opening_sdot_0(instance)* & + prm%sdot0* & (max(0.0_pReal,traction_n)/traction_crit - & - max(0.0_pReal,traction_n)/kinematics_slipplane_opening_critLoad(f,instance))**kinematics_slipplane_opening_N(instance) + max(0.0_pReal,traction_n)/prm%critLoad(i))**prm%n if (abs(udotn) > tol_math_check) then Ld = Ld + udotn*projection_n - dudotn_dt = udotn*kinematics_slipplane_opening_N(instance)/traction_n + dudotn_dt = udotn*prm%n/traction_n forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + & dudotn_dt*projection_n(k,l)*projection_n(m,n) endif - enddo enddo +end associate + end subroutine kinematics_slipplane_opening_LiAndItsTangent end module kinematics_slipplane_opening diff --git a/src/source_damage_anisoDuctile.f90 b/src/source_damage_anisoDuctile.f90 index dd5b95893..07f8e5e58 100644 --- a/src/source_damage_anisoDuctile.f90 +++ b/src/source_damage_anisoDuctile.f90 @@ -22,9 +22,6 @@ module source_damage_anisoDuctile source_damage_anisoDuctile_output !< name of each post result output - integer(pInt), dimension(:,:), allocatable, private :: & - source_damage_anisoDuctile_Nslip !< number of slip systems per family - enum, bind(c) enumerator :: undefined_ID, & damage_drivingforce_ID @@ -37,9 +34,9 @@ module source_damage_anisoDuctile N real(pReal), dimension(:), allocatable :: & critPlasticStrain - integer(pInt) :: & + integer :: & totalNslip - integer(pInt), dimension(:), allocatable :: & + integer, dimension(:), allocatable :: & Nslip integer(kind(undefined_ID)), allocatable, dimension(:) :: & outputID @@ -82,13 +79,10 @@ subroutine source_damage_anisoDuctile_init material_phase, & sourceState use config, only: & - config_phase, & - material_Nphase - use lattice, only: & - lattice_maxNslipFamily + config_phase + implicit none - integer(pInt) :: Ninstance,phase,instance,source,sourceOffset integer(pInt) :: NofMyPhase,p ,i @@ -110,9 +104,9 @@ subroutine source_damage_anisoDuctile_init if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - allocate(source_damage_anisoDuctile_offset(material_Nphase), source=0_pInt) - allocate(source_damage_anisoDuctile_instance(material_Nphase), source=0_pInt) - do phase = 1, material_Nphase + allocate(source_damage_anisoDuctile_offset(size(config_phase)), source=0_pInt) + allocate(source_damage_anisoDuctile_instance(size(config_phase)), source=0_pInt) + do phase = 1, size(config_phase) source_damage_anisoDuctile_instance(phase) = count(phase_source(:,1:phase) == source_damage_anisoDuctile_ID) do source = 1, phase_Nsources(phase) if (phase_source(source,phase) == source_damage_anisoDuctile_ID) & @@ -124,7 +118,6 @@ subroutine source_damage_anisoDuctile_init allocate(source_damage_anisoDuctile_output(maxval(phase_Noutput),Ninstance)) source_damage_anisoDuctile_output = '' - allocate(source_damage_anisoDuctile_Nslip(lattice_maxNslipFamily,Ninstance), source=0_pInt) allocate(param(Ninstance)) @@ -136,7 +129,7 @@ subroutine source_damage_anisoDuctile_init prm%aTol = config%getFloat('anisoductile_atol',defaultVal = 1.0e-3_pReal) prm%N = config%getFloat('anisoductile_ratesensitivity') - + prm%totalNslip = sum(prm%Nslip) ! sanity checks if (prm%aTol < 0.0_pReal) extmsg = trim(extmsg)//' anisoductile_atol' @@ -185,8 +178,6 @@ subroutine source_damage_anisoDuctile_init sourceState(phase)%p(sourceOffset)%sizePostResults = sum(source_damage_anisoDuctile_sizePostResult(:,instance)) sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol - source_damage_anisoDuctile_Nslip(1:size(param(instance)%Nslip),instance) = param(instance)%Nslip - enddo end subroutine source_damage_anisoDuctile_init @@ -202,8 +193,6 @@ subroutine source_damage_anisoDuctile_dotState(ipc, ip, el) material_homog, & damage, & damageMapping - use lattice, only: & - lattice_maxNslipFamily implicit none integer(pInt), intent(in) :: & @@ -216,7 +205,7 @@ subroutine source_damage_anisoDuctile_dotState(ipc, ip, el) sourceOffset, & homog, damageOffset, & instance, & - index, f, i + f, i phase = phaseAt(ipc,ip,el) constituent = phasememberAt(ipc,ip,el) @@ -225,17 +214,12 @@ subroutine source_damage_anisoDuctile_dotState(ipc, ip, el) homog = material_homog(ip,el) damageOffset = damageMapping(homog)%p(ip,el) - index = 1_pInt - sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = 0.0_pReal - do f = 1_pInt,lattice_maxNslipFamily - do i = 1_pInt,source_damage_anisoDuctile_Nslip(f,instance) ! process each (active) slip system in family + + do i = 1, param(instance)%totalNslip sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = & sourceState(phase)%p(sourceOffset)%dotState(1,constituent) + & - plasticState(phase)%slipRate(index,constituent)/ & - ((damage(homog)%p(damageOffset))**param(instance)%N)/param(instance)%critPlasticStrain(index) - - index = index + 1_pInt - enddo + plasticState(phase)%slipRate(i,constituent)/ & + ((damage(homog)%p(damageOffset))**param(instance)%N)/param(instance)%critPlasticStrain(i) enddo end subroutine source_damage_anisoDuctile_dotState