diff --git a/CMakeLists.txt b/CMakeLists.txt index ee376ef02..f742034e2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -280,7 +280,7 @@ if (CMAKE_Fortran_COMPILER_ID STREQUAL "Intel") # ... for uninitialized variables. set (DEBUG_FLAGS "${DEBUG_FLAGS} -ftrapuv") # ... initializes stack local variables to an unusual value to aid error detection - set (DEBUG_FLAGS "${DEBUG_FLAGS} -fpe-all0") + set (DEBUG_FLAGS "${DEBUG_FLAGS} -fpe-all=0") # ... capture all floating-point exceptions, sets -ftz automatically set (DEBUG_FLAGS "${DEBUG_FLAGS} -warn") diff --git a/src/material.f90 b/src/material.f90 index d54659acc..fc27b0cf4 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -283,6 +283,7 @@ module material public :: & material_init, & + material_allocatePlasticState, & ELASTICITY_hooke_ID ,& PLASTICITY_none_ID, & PLASTICITY_isotropic_ID, & @@ -1069,6 +1070,55 @@ subroutine material_parseTexture end subroutine material_parseTexture +!-------------------------------------------------------------------------------------------------- +!> @brief allocates the plastic state of a phase +!-------------------------------------------------------------------------------------------------- +subroutine material_allocatePlasticState(phase,NofMyPhase,sizeState,sizeDotState,sizeDeltaState,& + Nslip,Ntwin,Ntrans) + use numerics, only: & + numerics_integrator2 => numerics_integrator ! compatibility hack + + implicit none + integer(pInt), intent(in) :: & + phase, & + NofMyPhase, & + sizeState, & + sizeDotState, & + sizeDeltaState, & + Nslip, & + Ntwin, & + Ntrans + integer(pInt) :: numerics_integrator ! compatibility hack + numerics_integrator = numerics_integrator2(1) ! compatibility hack + + plasticState(phase)%sizeState = sizeState + plasticState(phase)%sizeDotState = sizeDotState + plasticState(phase)%sizeDeltaState = sizeDeltaState + plasticState(phase)%Nslip = Nslip + plasticState(phase)%Ntwin = Ntwin + plasticState(phase)%Ntrans= Ntrans + + allocate(plasticState(phase)%aTolState (sizeState), source=0.0_pReal) + allocate(plasticState(phase)%state0 (sizeState,NofMyPhase), source=0.0_pReal) + allocate(plasticState(phase)%partionedState0 (sizeState,NofMyPhase), source=0.0_pReal) + allocate(plasticState(phase)%subState0 (sizeState,NofMyPhase), source=0.0_pReal) + allocate(plasticState(phase)%state (sizeState,NofMyPhase), source=0.0_pReal) + + allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) + if (numerics_integrator == 1_pInt) then + allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) + allocate(plasticState(phase)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal) + endif + if (numerics_integrator == 4_pInt) & + allocate(plasticState(phase)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal) + if (numerics_integrator == 5_pInt) & + allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase), source=0.0_pReal) + + allocate(plasticState(phase)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) + +end subroutine material_allocatePlasticState + + !-------------------------------------------------------------------------------------------------- !> @brief populates the grains !> @details populates the grains by identifying active microstructure/homogenization pairs, diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 index 9d30cf184..57d48d109 100644 --- a/src/plastic_phenopowerlaw.f90 +++ b/src/plastic_phenopowerlaw.f90 @@ -32,10 +32,7 @@ module plastic_phenopowerlaw totalvolfrac_twin_ID end enum - type, private :: tParameters !< container type for internal constitutive parameters - integer(pInt) :: & - totalNslip, & - totalNtwin + type, private :: tParameters real(pReal) :: & gdot0_slip, & !< reference shear strain rate for slip gdot0_twin, & !< reference shear strain rate for twin @@ -50,45 +47,48 @@ module plastic_phenopowerlaw h0_TwinSlip, & !< reference hardening twin - slip h0_TwinTwin, & !< reference hardening twin - twin a_slip, & - aTolResistance, & ! default absolute tolerance 1 Pa - aTolShear, & ! default absolute tolerance 1e-6 - aTolTwinfrac ! default absolute tolerance 1e-6 - integer(pInt), dimension(:), allocatable :: & - Nslip, & !< active number of slip systems per family - Ntwin !< active number of twin systems per family - real(pReal), dimension(:), allocatable :: & + 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(:) :: & 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) !ToDo: Better name! gamma_twin_char !< characteristic shear for twins - real(pReal), dimension(:,:), allocatable :: & + 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), dimension(:,:,:), allocatable :: & + real(pReal), allocatable, dimension(:,:,:) :: & Schmid_slip, & Schmid_twin, & nonSchmid_pos, & nonSchmid_neg - integer(kind(undefined_ID)), dimension(:), allocatable :: & - outputID !< ID of each post result output - end type + integer(pInt) :: & + totalNslip, & !< total number of active slip system + totalNtwin !< total number of active twin systems + integer(pInt), 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 + end type !< container type for internal constitutive parameters type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) type, private :: tPhenopowerlawState - real(pReal), pointer, dimension(:,:) :: & + real(pReal), pointer, dimension(:) :: & + sumGamma, & ! ToDo: why not make a dependent state? + sumF ! ToDo: why not make a dependent state? + real(pReal), pointer, dimension(:,:) :: & xi_slip, & xi_twin, & gamma_slip, & gamma_twin, & whole - real(pReal), pointer, dimension(:) :: & - sumGamma, & - sumF end type type(tPhenopowerlawState), allocatable, dimension(:), private :: & @@ -115,6 +115,7 @@ subroutine plastic_phenopowerlaw_init compiler_options #endif use prec, only: & + pStringLen, & dEq0 use debug, only: & debug_level, & @@ -130,6 +131,7 @@ subroutine plastic_phenopowerlaw_init phase_plasticity, & phase_plasticityInstance, & phase_Noutput, & + material_allocatePlasticState, & PLASTICITY_PHENOPOWERLAW_LABEL, & PLASTICITY_PHENOPOWERLAW_ID, & material_phase, & @@ -138,20 +140,17 @@ subroutine plastic_phenopowerlaw_init MATERIAL_partPhase, & config_phase use lattice - use numerics,only: & - numerics_integrator implicit none - integer(pInt) :: & - maxNinstance, & - instance,p,j,k, o, i,& + Ninstance, & + p, i, & NipcMyPhase, outputSize, & - sizeState,sizeDotState, & + sizeState, sizeDotState, & startIndex, endIndex - integer(pInt), dimension(0), parameter :: emptyIntArray = [integer(pInt)::] - real(pReal), dimension(0), parameter :: emptyRealArray = [real(pReal)::] + integer(pInt), dimension(0), parameter :: emptyIntArray = [integer(pInt)::] + real(pReal), dimension(0), parameter :: emptyRealArray = [real(pReal)::] character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] type(tParameters) :: & @@ -163,152 +162,151 @@ subroutine plastic_phenopowerlaw_init integer(kind(undefined_ID)) :: & outputID !< ID of each post result output - character(len=512) :: & - extmsg = '', & - structure = '' - character(len=65536), dimension(:), allocatable :: outputs + character(len=pStringLen) :: & + structure = '',& + extmsg = '' + character(len=65536), dimension(:), allocatable :: & + outputs write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_PHENOPOWERLAW_label//' init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - maxNinstance = int(count(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID),pInt) + Ninstance = int(count(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID),pInt) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & - write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance + write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - allocate(plastic_phenopowerlaw_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt) - allocate(plastic_phenopowerlaw_output(maxval(phase_Noutput),maxNinstance)) + allocate(plastic_phenopowerlaw_sizePostResult(maxval(phase_Noutput),Ninstance),source=0_pInt) + allocate(plastic_phenopowerlaw_output(maxval(phase_Noutput),Ninstance)) plastic_phenopowerlaw_output = '' - allocate(param(maxNinstance)) ! one container of parameters per instance - allocate(state(maxNinstance)) - allocate(dotState(maxNinstance)) + allocate(param(Ninstance)) + allocate(state(Ninstance)) + allocate(dotState(Ninstance)) do p = 1_pInt, size(phase_plasticityInstance) if (phase_plasticity(p) /= PLASTICITY_PHENOPOWERLAW_ID) cycle - instance = phase_plasticityInstance(p) - associate(prm => param(instance),stt => state(instance),dot => dotState(instance)) - extmsg = '' + associate(prm => param(phase_plasticityInstance(p)), & + dot => dotState(phase_plasticityInstance(p)), & + stt => state(phase_plasticityInstance(p))) - structure = config_phase(p)%getString('lattice_structure') + structure = config_phase(p)%getString('lattice_structure') +!-------------------------------------------------------------------------------------------------- +! optional parameters that need to be defined + prm%twinB = config_phase(p)%getFloat('twin_b',defaultVal=1.0_pReal) + prm%twinC = config_phase(p)%getFloat('twin_c',defaultVal=0.0_pReal) + prm%twinD = config_phase(p)%getFloat('twin_d',defaultVal=0.0_pReal) + prm%twinE = config_phase(p)%getFloat('twin_e',defaultVal=0.0_pReal) + + prm%aTolResistance = config_phase(p)%getFloat('atol_resistance',defaultVal=1.0_pReal) + prm%aTolShear = config_phase(p)%getFloat('atol_shear', defaultVal=1.0e-6_pReal) + prm%aTolTwinfrac = config_phase(p)%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_phase(p)%getInts('nslip',defaultVal=emptyIntArray) prm%totalNslip = sum(prm%Nslip) - if (size(prm%Nslip) > count(lattice_NslipSystem(:,p) > 0_pInt)) & - call IO_error(150_pInt,ext_msg='Nslip') - if (any(lattice_NslipSystem(1:size(prm%Nslip),p)-prm%Nslip < 0_pInt)) & - call IO_error(150_pInt,ext_msg='Nslip') - slipActive: if (prm%totalNslip > 0_pInt) then - prm%Schmid_slip = lattice_SchmidMatrix_slip(prm%Nslip,structure(1:3),& config_phase(p)%getFloat('c/a',defaultVal=0.0_pReal)) - ! reading in slip related parameters + if(structure=='bcc') then + prm%nonSchmidCoeff = config_phase(p)%getFloats('nonschmid_coefficients',& + defaultVal = emptyRealArray) + prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt) + prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1_pInt) + else + prm%nonSchmid_pos = prm%Schmid_slip + prm%nonSchmid_neg = prm%Schmid_slip + endif + prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip, & + config_phase(p)%getFloats('interaction_slipslip'), & + structure(1:3)) + prm%xi_slip_0 = config_phase(p)%getFloats('tau0_slip', requiredShape=shape(prm%Nslip)) prm%xi_slip_sat = config_phase(p)%getFloats('tausat_slip', requiredShape=shape(prm%Nslip)) - prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip,config_phase(p)%getFloats('interaction_slipslip'), & - structure(1:3)) prm%H_int = config_phase(p)%getFloats('h_int', requiredShape=shape(prm%Nslip), & defaultVal=[(0.0_pReal,i=1_pInt,size(prm%Nslip))]) - prm%nonSchmidCoeff = config_phase(p)%getFloats('nonschmid_coefficients',& - defaultVal = emptyRealArray ) - if(structure=='bcc') then - prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt) - prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1_pInt) - else - prm%nonSchmid_pos = prm%Schmid_slip - prm%nonSchmid_neg = prm%Schmid_slip - endif + prm%gdot0_slip = config_phase(p)%getFloat('gdot0_slip') prm%n_slip = config_phase(p)%getFloat('n_slip') prm%a_slip = config_phase(p)%getFloat('a_slip') prm%h0_SlipSlip = config_phase(p)%getFloat('h0_slipslip') - ! sanity checks for slip related parameters - if (any(prm%xi_slip_0 < 0.0_pReal .and. prm%Nslip > 0_pInt)) & - extmsg = trim(extmsg)//"xi_slip_0 " - if (any(prm%xi_slip_sat < prm%xi_slip_0 .and. prm%Nslip > 0_pInt)) & - extmsg = trim(extmsg)//"xi_slip_sat " - - if (prm%gdot0_slip <= 0.0_pReal) extmsg = trim(extmsg)//"gdot0_slip " - if (dEq0(prm%a_slip)) extmsg = trim(extmsg)//"a_slip " ! ToDo: negative values ok? - if (dEq0(prm%n_slip)) extmsg = trim(extmsg)//"n_slip " ! ToDo: negative values ok? - - ! expand slip related parameters from system => family - prm%xi_slip_0 = math_expand(prm%xi_slip_0,prm%Nslip) + ! 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) + 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 (dEq0(prm%a_slip)) extmsg = trim(extmsg)//'a_slip ' ! ToDo: negative values ok? + if (dEq0(prm%n_slip)) extmsg = trim(extmsg)//'n_slip ' ! ToDo: negative values ok? + if (any(prm%xi_slip_0 <= 0.0_pReal)) extmsg = trim(extmsg)//'xi_slip_0 ' + if (any(prm%xi_slip_sat < prm%xi_slip_0)) extmsg = trim(extmsg)//'xi_slip_sat ' else slipActive allocate(prm%interaction_SlipSlip(0,0)) allocate(prm%xi_slip_0(0)) endif slipActive +!-------------------------------------------------------------------------------------------------- +! twin related parameters prm%Ntwin = config_phase(p)%getInts('ntwin', defaultVal=emptyIntArray) prm%totalNtwin = sum(prm%Ntwin) - if (size(prm%Ntwin) > count(lattice_NtwinSystem(:,p) > 0_pInt)) & - call IO_error(150_pInt,ext_msg='Ntwin') - if (any(lattice_NtwinSystem(1:size(prm%Ntwin),p)-prm%Ntwin < 0_pInt)) & - call IO_error(150_pInt,ext_msg='Ntwin') - twinActive: if (prm%totalNtwin > 0_pInt) then prm%Schmid_twin = lattice_SchmidMatrix_twin(prm%Ntwin,structure(1:3),& config_phase(p)%getFloat('c/a',defaultVal=0.0_pReal)) - ! reading in twin related parameters + prm%interaction_TwinTwin = lattice_interaction_TwinTwin(prm%Ntwin,& + config_phase(p)%getFloats('interaction_twintwin'), & + structure(1:3)) + prm%gamma_twin_char = lattice_characteristicShear_twin(prm%Ntwin,structure(1:3),& + config_phase(p)%getFloat('c/a')) + prm%xi_twin_0 = config_phase(p)%getFloats('tau0_twin',requiredShape=shape(prm%Ntwin)) - prm%interaction_TwinTwin = lattice_interaction_TwinTwin(prm%Ntwin,config_phase(p)%getFloats('interaction_twintwin'), & - structure(1:3)) - prm%gdot0_twin = config_phase(p)%getFloat('gdot0_twin') - prm%n_twin = config_phase(p)%getFloat('n_twin') - prm%spr = config_phase(p)%getFloat('s_pr') - prm%h0_TwinTwin = config_phase(p)%getFloat('h0_twintwin') + prm%gdot0_twin = config_phase(p)%getFloat('gdot0_twin') + prm%n_twin = config_phase(p)%getFloat('n_twin') + prm%spr = config_phase(p)%getFloat('s_pr') + prm%h0_TwinTwin = config_phase(p)%getFloat('h0_twintwin') - ! sanity checks for twin related parameters - if (any(prm%xi_twin_0 < 0.0_pReal .and. prm%Ntwin > 0_pInt)) & - extmsg = trim(extmsg)//"xi_twin_0 " - if (prm%gdot0_twin <= 0.0_pReal) extmsg = trim(extmsg)//"gdot0_twin " - if (dEq0(prm%n_twin)) extmsg = trim(extmsg)//"n_twin " ! ToDo: negative values ok? + ! expand: family => system + prm%xi_twin_0 = math_expand(prm%xi_twin_0, prm%Ntwin) - ! expand slip related parameters from system => family - 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 (dEq0(prm%n_twin)) extmsg = trim(extmsg)//'n_twin ' ! ToDo: negative values ok? else twinActive allocate(prm%interaction_TwinTwin(0,0)) allocate(prm%xi_twin_0(0)) endif twinActive - prm%gamma_twin_char = lattice_characteristicShear_twin(prm%Ntwin,structure(1:3),& - config_phase(p)%getFloat('c/a',defaultVal=0.0_pReal)) - +!-------------------------------------------------------------------------------------------------- +! slip-twin related parameters slipAndTwinActive: if (prm%totalNslip > 0_pInt .and. prm%totalNtwin > 0_pInt) then prm%interaction_SlipTwin = lattice_interaction_SlipTwin(prm%Nslip,prm%Ntwin,& - config_phase(p)%getFloats('interaction_sliptwin'), & - structure(1:3)) + config_phase(p)%getFloats('interaction_sliptwin'), & + structure(1:3)) prm%interaction_TwinSlip = lattice_interaction_TwinSlip(prm%Ntwin,prm%Nslip,& - config_phase(p)%getFloats('interaction_twinslip'), & - structure(1:3)) + config_phase(p)%getFloats('interaction_twinslip'), & + structure(1:3)) else slipAndTwinActive - allocate(prm%interaction_SlipTwin(prm%totalNslip,prm%TotalNtwin)) ! at least one dimension 0 - allocate(prm%interaction_TwinSlip(prm%totalNtwin,prm%TotalNslip)) ! at least one dimension 0 + allocate(prm%interaction_SlipTwin(prm%totalNslip,prm%TotalNtwin)) ! at least one dimension is 0 + allocate(prm%interaction_TwinSlip(prm%totalNtwin,prm%TotalNslip)) ! at least one dimension is 0 prm%h0_TwinSlip = 0.0_pReal endif slipAndTwinActive - ! optional parameters that should be defined - prm%twinB = config_phase(p)%getFloat('twin_b',defaultVal=1.0_pReal) - prm%twinC = config_phase(p)%getFloat('twin_c',defaultVal=0.0_pReal) - prm%twinD = config_phase(p)%getFloat('twin_d',defaultVal=0.0_pReal) - prm%twinE = config_phase(p)%getFloat('twin_e',defaultVal=0.0_pReal) - - prm%aTolResistance = config_phase(p)%getFloat('atol_resistance',defaultVal=1.0_pReal) - prm%aTolShear = config_phase(p)%getFloat('atol_shear', defaultVal=1.0e-6_pReal) - prm%aTolTwinfrac = config_phase(p)%getFloat('atol_twinfrac', defaultVal=1.0e-6_pReal) - - 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 " - - if (extmsg /= '') call IO_error(211_pInt,ip=instance,& - ext_msg=trim(extmsg)//'('//PLASTICITY_PHENOPOWERLAW_label//')') +!-------------------------------------------------------------------------------------------------- +! exit if any parameter is out of range + if (extmsg /= '') & + call IO_error(211_pInt,ext_msg=trim(extmsg)//'('//PLASTICITY_PHENOPOWERLAW_label//')') +!-------------------------------------------------------------------------------------------------- +! output pararameters outputs = config_phase(p)%getStrings('(output)',defaultVal=emptyStringArray) allocate(prm%outputID(0)) do i=1_pInt, size(outputs) @@ -349,8 +347,8 @@ subroutine plastic_phenopowerlaw_init end select if (outputID /= undefined_ID) then - plastic_phenopowerlaw_output(i,instance) = outputs(i) - plastic_phenopowerlaw_sizePostResult(i,instance) = outputSize + plastic_phenopowerlaw_output(i,phase_plasticityInstance(p)) = outputs(i) + plastic_phenopowerlaw_sizePostResult(i,phase_plasticityInstance(p)) = outputSize prm%outputID = [prm%outputID , outputID] endif @@ -361,32 +359,12 @@ subroutine plastic_phenopowerlaw_init NipcMyPhase = count(material_phase == p) ! number of IPCs containing my phase sizeState = size(['tau_slip ','gamma_slip']) * prm%TotalNslip & + size(['tau_twin ','gamma_twin']) * prm%TotalNtwin & - + size(['sum(gamma)','sum(f) ']) - -!-------------------------------------------------------------------------------------------------- -! ToDo: This could be done by a function (in constitutive?) + + size(['sum(gamma)','sum(f) ']) ! ToDo: only needed if either twin or slip active! sizeDotState = sizeState - plasticState(p)%sizeState = sizeState - plasticState(p)%sizeDotState = sizeDotState - plasticState(p)%sizePostResults = sum(plastic_phenopowerlaw_sizePostResult(:,instance)) - plasticState(p)%nSlip = prm%totalNslip - plasticState(p)%nTwin = prm%totalNtwin - allocate(plasticState(p)%aTolState ( sizeState), source=0.0_pReal) - allocate(plasticState(p)%state0 ( sizeState,NipcMyPhase), source=0.0_pReal) - allocate(plasticState(p)%partionedState0 ( sizeState,NipcMyPhase), source=0.0_pReal) - allocate(plasticState(p)%subState0 ( sizeState,NipcMyPhase), source=0.0_pReal) - allocate(plasticState(p)%state ( sizeState,NipcMyPhase), source=0.0_pReal) - allocate(plasticState(p)%dotState (sizeDotState,NipcMyPhase), source=0.0_pReal) - allocate(plasticState(p)%deltaState (0_pInt,NipcMyPhase), source=0.0_pReal) - if (any(numerics_integrator == 1_pInt)) then - allocate(plasticState(p)%previousDotState (sizeDotState,NipcMyPhase),source=0.0_pReal) - allocate(plasticState(p)%previousDotState2(sizeDotState,NipcMyPhase),source=0.0_pReal) - endif - if (any(numerics_integrator == 4_pInt)) & - allocate(plasticState(p)%RK4dotState (sizeDotState,NipcMyPhase), source=0.0_pReal) - if (any(numerics_integrator == 5_pInt)) & - allocate(plasticState(p)%RKCK45dotState (6,sizeDotState,NipcMyPhase), source=0.0_pReal) + call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0_pInt, & + prm%totalNslip,prm%totalNtwin,0_pInt) + plasticState(p)%sizePostResults = sum(plastic_phenopowerlaw_sizePostResult(:,phase_plasticityInstance(p))) !-------------------------------------------------------------------------------------------------- @@ -432,7 +410,7 @@ subroutine plastic_phenopowerlaw_init dot%gamma_twin => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear - plasticState(p)%state0 = plasticState(p)%state + plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally dot%whole => plasticState(p)%dotState end associate @@ -544,25 +522,21 @@ subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) dot%gamma_slip(:,of) = abs(gdot_slip_pos+gdot_slip_neg) dot%sumGamma(of) = sum(dot%gamma_slip(:,of)) call kinetics_twin(prm,stt,of,Mp,dot%gamma_twin(:,of)) - if (stt%sumF(of) < 0.98_pReal) dot%sumF(of) = sum(dot%gamma_twin(:,of)/prm%gamma_twin_char) + if (prm%totalNtwin > 0_pInt) dot%sumF(of) = merge(sum(dot%gamma_twin(:,of)/prm%gamma_twin_char), & + 0.0_pReal, & + stt%sumF(of) < 0.98_pReal) !-------------------------------------------------------------------------------------------------- ! hardening hardeningSlip: do i = 1_pInt, prm%totalNslip - dot%xi_slip(i,of) = & - c_SlipSlip * left_SlipSlip(i) & - * dot_product(prm%interaction_SlipSlip(i,:),right_SlipSlip*dot%gamma_slip(:,of)) & - + & - dot_product(prm%interaction_SlipTwin(i,:),dot%gamma_twin(:,of)) + dot%xi_slip(i,of) = dot_product(prm%interaction_SlipSlip(i,:),right_SlipSlip*dot%gamma_slip(:,of)) & + * c_SlipSlip * left_SlipSlip(i) & + + dot_product(prm%interaction_SlipTwin(i,:),dot%gamma_twin(:,of)) enddo hardeningSlip hardeningTwin: do i = 1_pInt, prm%totalNtwin - dot%xi_twin(i,of) = & - c_TwinSlip & - * dot_product(prm%interaction_TwinSlip(i,:),dot%gamma_slip(:,of)) & - + & - c_TwinTwin & - * dot_product(prm%interaction_TwinTwin(i,:),dot%gamma_twin(:,of)) + dot%xi_twin(i,of) = c_TwinSlip * dot_product(prm%interaction_TwinSlip(i,:),dot%gamma_slip(:,of)) & + + c_TwinTwin * dot_product(prm%interaction_TwinTwin(i,:),dot%gamma_twin(:,of)) enddo hardeningTwin end associate @@ -572,7 +546,7 @@ end subroutine plastic_phenopowerlaw_dotState !-------------------------------------------------------------------------------------------------- !> @brief calculates shear rates on slip systems and derivatives with respect to resolved stress -!> @details: Shear rates are calculated only optionally. NOTE: Agains the common convention, the +!> @details: Shear rates 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_slip(prm,stt,of,Mp,gdot_slip_pos,gdot_slip_neg, & @@ -602,26 +576,39 @@ pure subroutine kinetics_slip(prm,stt,of,Mp,gdot_slip_pos,gdot_slip_neg, & tau_slip_pos, & tau_slip_neg integer(pInt) :: i + logical :: nonSchmidActive + + nonSchmidActive = size(prm%nonSchmidCoeff) > 0_pInt do i = 1_pInt, prm%totalNslip - tau_slip_pos(i) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i)) - tau_slip_neg(i) = math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i)) + 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 - gdot_slip_pos = 0.5_pReal*prm%gdot0_slip & - * sign(abs(tau_slip_pos/stt%xi_slip(:,of))**prm%n_slip, tau_slip_pos) - gdot_slip_neg = 0.5_pReal*prm%gdot0_slip & - * sign(abs(tau_slip_neg/stt%xi_slip(:,of))**prm%n_slip, tau_slip_neg) + 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 = 0.5_pReal*prm%gdot0_slip & + * 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(tau_slip_pos)) + where(dNeq0(gdot_slip_pos)) dgdot_dtau_slip_pos = gdot_slip_pos*prm%n_slip/tau_slip_pos else where dgdot_dtau_slip_pos = 0.0_pReal end where endif if (present(dgdot_dtau_slip_neg)) then - where(dNeq0(tau_slip_neg)) + where(dNeq0(gdot_slip_neg)) dgdot_dtau_slip_neg = gdot_slip_neg*prm%n_slip/tau_slip_neg else where dgdot_dtau_slip_neg = 0.0_pReal @@ -633,7 +620,7 @@ end subroutine kinetics_slip !-------------------------------------------------------------------------------------------------- !> @brief calculates shear rates on twin systems and derivatives with respect to resolved stress -!> @details: Shear rates are calculated only optionally. NOTE: Agains the common convention, the +!> @details: Shear rates 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(prm,stt,of,Mp,gdot_twin,dgdot_dtau_twin) @@ -663,11 +650,15 @@ pure subroutine kinetics_twin(prm,stt,of,Mp,gdot_twin,dgdot_dtau_twin) do i = 1_pInt, prm%totalNtwin tau_twin(i) = math_mul33xx33(Mp,prm%Schmid_twin(1:3,1:3,i)) enddo - gdot_twin = merge((1.0_pReal-stt%sumF(of))*prm%gdot0_twin*(abs(tau_twin)/stt%xi_twin(:,of))**prm%n_twin, & - 0.0_pReal, tau_twin>0.0_pReal) + + where(tau_twin > 0.0_pReal) + gdot_twin = (1.0_pReal-stt%sumF(of))*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(tau_twin)) + where(dNeq0(gdot_twin)) dgdot_dtau_twin = gdot_twin*prm%n_twin/tau_twin else where dgdot_dtau_twin = 0.0_pReal @@ -681,14 +672,8 @@ end subroutine kinetics_twin !> @brief return array of constitutive results !-------------------------------------------------------------------------------------------------- function plastic_phenopowerlaw_postResults(Mp,instance,of) result(postResults) - use material, only: & - material_phase, & - plasticState, & - phasememberAt, & - phase_plasticityInstance use math, only: & - math_mul33xx33, & - math_Mandel6to33 + math_mul33xx33 implicit none real(pReal), dimension(3,3), intent(in) :: & @@ -701,7 +686,7 @@ function plastic_phenopowerlaw_postResults(Mp,instance,of) result(postResults) postResults integer(pInt) :: & - o,c,i,j + o,c,i real(pReal), dimension(param(instance)%totalNslip) :: & gdot_slip_pos,gdot_slip_neg