diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index eb1448028..6b32e39ec 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -74,6 +74,7 @@ add_library (PLASTIC OBJECT "plastic_disloUCLA.f90" "plastic_isotropic.f90" "plastic_phenopowerlaw.f90" + "plastic_MITdiag.f90" "plastic_titanmod.f90" "plastic_nonlocal.f90" "plastic_none.f90" diff --git a/src/constitutive.f90 b/src/constitutive.f90 index de8f61c2a..99594a846 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -70,6 +70,7 @@ subroutine constitutive_init() PLASTICITY_none_ID, & PLASTICITY_isotropic_ID, & PLASTICITY_phenopowerlaw_ID, & + PLASTICITY_mitdiag_ID ,& PLASTICITY_phenoplus_ID, & PLASTICITY_dislotwin_ID, & PLASTICITY_disloucla_ID, & @@ -93,6 +94,7 @@ subroutine constitutive_init() PLASTICITY_NONE_label, & PLASTICITY_ISOTROPIC_label, & PLASTICITY_PHENOPOWERLAW_label, & + PLASTICITY_MITDIAG_label, & PLASTICITY_PHENOPLUS_label, & PLASTICITY_DISLOTWIN_label, & PLASTICITY_DISLOUCLA_label, & @@ -113,6 +115,7 @@ subroutine constitutive_init() use plastic_none use plastic_isotropic use plastic_phenopowerlaw + use plastic_mitdiag use plastic_phenoplus use plastic_dislotwin use plastic_disloucla @@ -158,6 +161,7 @@ subroutine constitutive_init() if (any(phase_plasticity == PLASTICITY_NONE_ID)) call plastic_none_init if (any(phase_plasticity == PLASTICITY_ISOTROPIC_ID)) call plastic_isotropic_init(FILEUNIT) if (any(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)) call plastic_phenopowerlaw_init(FILEUNIT) + if (any(phase_plasticity == PLASTICITY_MITDIAG_ID)) call plastic_mitdiag_init(FILEUNIT) if (any(phase_plasticity == PLASTICITY_PHENOPLUS_ID)) call plastic_phenoplus_init(FILEUNIT) if (any(phase_plasticity == PLASTICITY_DISLOTWIN_ID)) call plastic_dislotwin_init(FILEUNIT) if (any(phase_plasticity == PLASTICITY_DISLOUCLA_ID)) call plastic_disloucla_init(FILEUNIT) @@ -207,7 +211,7 @@ subroutine constitutive_init() outputName = PLASTICITY_NONE_label thisNoutput => null() thisOutput => null() - thisSize => null() + thisSize => null() case (PLASTICITY_ISOTROPIC_ID) plasticityType outputName = PLASTICITY_ISOTROPIC_label thisNoutput => plastic_isotropic_Noutput @@ -218,6 +222,11 @@ subroutine constitutive_init() thisNoutput => plastic_phenopowerlaw_Noutput thisOutput => plastic_phenopowerlaw_output thisSize => plastic_phenopowerlaw_sizePostResult + case (PLASTICITY_MITDIAG_ID) plasticityType + outputName = PLASTICITY_MITDIAG_label + thisNoutput => plastic_mitdiag_Noutput + thisOutput => plastic_mitdiag_output + thisSize => plastic_mitdiag_sizePostResult case (PLASTICITY_PHENOPLUS_ID) plasticityType outputName = PLASTICITY_PHENOPLUS_label thisNoutput => plastic_phenoplus_Noutput @@ -501,6 +510,7 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v PLASTICITY_NONE_ID, & PLASTICITY_ISOTROPIC_ID, & PLASTICITY_PHENOPOWERLAW_ID, & + PLASTICITY_MITDIAG_ID, & PLASTICITY_PHENOPLUS_ID, & PLASTICITY_DISLOTWIN_ID, & PLASTICITY_DISLOUCLA_ID, & @@ -510,6 +520,8 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v plastic_isotropic_LpAndItsTangent use plastic_phenopowerlaw, only: & plastic_phenopowerlaw_LpAndItsTangent + use plastic_mitdiag, only: & + plastic_mitdiag_LpAndItsTangent use plastic_phenoplus, only: & plastic_phenoplus_LpAndItsTangent use plastic_dislotwin, only: & @@ -560,6 +572,8 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v call plastic_isotropic_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v,ipc,ip,el) case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType call plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v,ipc,ip,el) + case (PLASTICITY_MITDIAG_ID) plasticityType + call plastic_mitdiag_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v,ipc,ip,el) case (PLASTICITY_PHENOPLUS_ID) plasticityType call plastic_phenoplus_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v,ipc,ip,el) case (PLASTICITY_NONLOCAL_ID) plasticityType @@ -884,6 +898,7 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra PLASTICITY_none_ID, & PLASTICITY_isotropic_ID, & PLASTICITY_phenopowerlaw_ID, & + PLASTICITY_mitdiag_ID, & PLASTICITY_phenoplus_ID, & PLASTICITY_dislotwin_ID, & PLASTICITY_disloucla_ID, & @@ -897,6 +912,8 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra plastic_isotropic_dotState use plastic_phenopowerlaw, only: & plastic_phenopowerlaw_dotState + use plastic_mitdiag, only: & + plastic_mitdiag_dotState use plastic_phenoplus, only: & plastic_phenoplus_dotState use plastic_dislotwin, only: & @@ -950,6 +967,8 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra call plastic_isotropic_dotState (Tstar_v,ipc,ip,el) case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType call plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el) + case (PLASTICITY_MITDIAG_ID) plasticityType + call plastic_mitdiag_dotState(Tstar_v,ipc,ip,el) case (PLASTICITY_PHENOPLUS_ID) plasticityType call plastic_phenoplus_dotState (Tstar_v,ipc,ip,el) case (PLASTICITY_DISLOTWIN_ID) plasticityType @@ -1093,6 +1112,7 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el) PLASTICITY_NONE_ID, & PLASTICITY_ISOTROPIC_ID, & PLASTICITY_PHENOPOWERLAW_ID, & + PLASTICITY_MITDIAG_ID, & PLASTICITY_PHENOPLUS_ID, & PLASTICITY_DISLOTWIN_ID, & PLASTICITY_DISLOUCLA_ID, & @@ -1108,6 +1128,8 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el) plastic_phenopowerlaw_postResults use plastic_phenoplus, only: & plastic_phenoplus_postResults + use plastic_mitdiag, only: & + plastic_mitdiag_postResults use plastic_dislotwin, only: & plastic_dislotwin_postResults use plastic_disloucla, only: & @@ -1160,6 +1182,9 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el) case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType constitutive_postResults(startPos:endPos) = & plastic_phenopowerlaw_postResults(Tstar_v,ipc,ip,el) + case (PLASTICITY_MITDIAG_ID) plasticityType + constitutive_postResults(startPos:endPos) = & + plastic_mitdiag_postResults(Tstar_v,ipc,ip,el) case (PLASTICITY_PHENOPLUS_ID) plasticityType constitutive_postResults(startPos:endPos) = & plastic_phenoplus_postResults(Tstar_v,ipc,ip,el) diff --git a/src/material.f90 b/src/material.f90 index a77c4871a..1acff4f93 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -25,6 +25,7 @@ module material PLASTICITY_none_label = 'none', & PLASTICITY_isotropic_label = 'isotropic', & PLASTICITY_phenopowerlaw_label = 'phenopowerlaw', & + PLASTICITY_mitdiag_label = 'mitdiag', & PLASTICITY_phenoplus_label = 'phenoplus', & PLASTICITY_dislotwin_label = 'dislotwin', & PLASTICITY_disloucla_label = 'disloucla', & @@ -73,6 +74,7 @@ module material enumerator :: PLASTICITY_undefined_ID, & PLASTICITY_none_ID, & PLASTICITY_isotropic_ID, & + PLASTICITY_mitdiag_ID, & PLASTICITY_phenopowerlaw_ID, & PLASTICITY_phenoplus_ID, & PLASTICITY_dislotwin_ID, & @@ -312,6 +314,7 @@ module material PLASTICITY_none_ID, & PLASTICITY_isotropic_ID, & PLASTICITY_phenopowerlaw_ID, & + PLASTICITY_mitdiag_ID, & PLASTICITY_phenoplus_ID, & PLASTICITY_dislotwin_ID, & PLASTICITY_disloucla_ID, & @@ -671,7 +674,7 @@ subroutine material_parseHomogenization(fileUnit,myPart) case ('porosity') select case (IO_lc(IO_stringValue(line,chunkPos,2_pInt))) - case(POROSITY_NONE_label) + case(POROSITY_none_label) porosity_type(section) = POROSITY_none_ID case(POROSITY_phasefield_label) porosity_type(section) = POROSITY_phasefield_ID @@ -729,8 +732,6 @@ end subroutine material_parseHomogenization !> @brief parses the microstructure part in the material configuration file !-------------------------------------------------------------------------------------------------- subroutine material_parseMicrostructure(fileUnit,myPart) - use prec, only: & - dNeq use IO use mesh, only: & mesh_element, & @@ -740,6 +741,7 @@ subroutine material_parseMicrostructure(fileUnit,myPart) character(len=*), intent(in) :: myPart integer(pInt), intent(in) :: fileUnit + integer(pInt), allocatable, dimension(:) :: chunkPos integer(pInt) :: Nsections, section, constituent, e, i character(len=65536) :: & @@ -759,7 +761,7 @@ subroutine material_parseMicrostructure(fileUnit,myPart) allocate(microstructure_elemhomo(Nsections), source=.false.) if(any(mesh_element(4,1:mesh_NcpElems) > Nsections)) & - call IO_error(155_pInt,ext_msg='More microstructures in geometry than sections in material.config') + call IO_error(155_pInt,ext_msg='Microstructure in geometry > Sections in material.config') forall (e = 1_pInt:mesh_NcpElems) microstructure_active(mesh_element(4,e)) = .true. ! current microstructure used in model? Elementwise view, maximum N operations for N elements @@ -802,7 +804,7 @@ subroutine material_parseMicrostructure(fileUnit,myPart) microstructure_crystallite(section) = IO_intValue(line,chunkPos,2_pInt) case ('(constituent)') constituent = constituent + 1_pInt - do i = 2_pInt,6_pInt,2_pInt + do i=2_pInt,6_pInt,2_pInt tag = IO_lc(IO_stringValue(line,chunkPos,i)) select case (tag) case('phase') @@ -817,12 +819,6 @@ subroutine material_parseMicrostructure(fileUnit,myPart) endif enddo - !sanity check -do section = 1_pInt, Nsections - if (dNeq(sum(microstructure_fraction(:,section)),1.0_pReal)) & - call IO_error(153_pInt,ext_msg=microstructure_name(section)) -enddo - end subroutine material_parseMicrostructure @@ -979,17 +975,19 @@ subroutine material_parsePhase(fileUnit,myPart) end select case ('plasticity') select case (IO_lc(IO_stringValue(line,chunkPos,2_pInt))) - case (PLASTICITY_NONE_label) - phase_plasticity(section) = PLASTICITY_NONE_ID + case (PLASTICITY_none_label) + phase_plasticity(section) = PLASTICITY_none_ID case (PLASTICITY_ISOTROPIC_label) phase_plasticity(section) = PLASTICITY_ISOTROPIC_ID + case (PLASTICITY_mitdiag_label) + phase_plasticity(section) = PLASTICITY_MITDIAG_ID case (PLASTICITY_PHENOPOWERLAW_label) phase_plasticity(section) = PLASTICITY_PHENOPOWERLAW_ID case (PLASTICITY_PHENOPLUS_label) phase_plasticity(section) = PLASTICITY_PHENOPLUS_ID case (PLASTICITY_DISLOTWIN_label) phase_plasticity(section) = PLASTICITY_DISLOTWIN_ID - case (PLASTICITY_DISLOUCLA_label) + case (PLASTICITY_disloUCLA_label) phase_plasticity(section) = PLASTICITY_DISLOUCLA_ID case (PLASTICITY_TITANMOD_label) phase_plasticity(section) = PLASTICITY_TITANMOD_ID diff --git a/src/plastic_MITdiag.f90 b/src/plastic_MITdiag.f90 new file mode 100644 index 000000000..28d9bc3d7 --- /dev/null +++ b/src/plastic_MITdiag.f90 @@ -0,0 +1,716 @@ +!-------------------------------------------------------------------------------------------------- +!> @author Tias Maiti, Michigan State University +!> @author Philip Eisenlohr, Michigan State University +!> @brief material subroutine implementing the diagonal hardening concept outlined by +!> Raul Radovitzky in the paper +!-------------------------------------------------------------------------------------------------- + +module plastic_MITdiag + use prec, only: & + pReal,& + pInt + use lattice, only: & + lattice_maxNinteraction, & + lattice_maxNslipFamily + + implicit none + private + integer(pInt), dimension(:), allocatable, public, protected :: & + plastic_MITdiag_sizePostResults !< cumulative size of post results + + integer(pInt), dimension(:,:), allocatable, target, public :: & + plastic_MITdiag_sizePostResult !< size of each post result output + + character(len=64), dimension(:,:), allocatable, target, public :: & + plastic_MITdiag_output !< name of each post result output + + integer(pInt), dimension(:), allocatable, target, public :: & + plastic_MITdiag_Noutput !< number of outputs per instance + + integer(pInt), dimension(:), allocatable, public, protected :: & + plastic_MITdiag_totalNslip !< no. of slip system used in simulation + + real(pReal), dimension(:,:,:), allocatable, private :: & + plastic_MITdiag_hardeningMatrix + + enum, bind(c) + enumerator :: undefined_ID, & + resistance_ID, & + accumulatedshear_ID, & + shearrate_ID, & + resolvedstress_ID + end enum + + type, private :: tParameters !< container type for internal constitutive parameters + integer(kind(undefined_ID)), allocatable, dimension(:) :: & + outputID + real(pReal), dimension(lattice_maxNinteraction) :: & + interaction + real(pReal), dimension(:,:), allocatable :: & + hardeningMatrix + integer(pInt), dimension(lattice_maxNslipFamily) :: & + n_slip = 0_pInt + real(pReal), dimension(lattice_maxNslipFamily) :: & + b, & + g0, & + rho0, & + rhosat, & + gammasat + real(pReal) :: & + aTolResistance = 1.0_pReal, & + aTolShear = 1.0e-6_pReal, & + gdot0, & + n + end type + + type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) + + type, private :: tMITdiagState !< internal state aliases + real(pReal), pointer, dimension(:,:) :: & ! scalars along NipcMyInstance + resistance, & + accumulatedShear + end type + type, private :: tMITdiagAbsTol !< internal alias for abs tolerance in state + real(pReal), pointer :: & ! scalars along NipcMyInstance + resistance, & + accumulatedShear + end type + type(tMITdiagState), allocatable, dimension(:), private :: & !< state aliases per instance + state, & + state0, & + dotState + type(tMITdiagAbsTol), allocatable, dimension(:), private :: & !< state aliases per instance + stateAbsTol + + public :: & + plastic_MITdiag_init, & + plastic_MITdiag_LpAndItsTangent, & + plastic_MITdiag_dotState, & + plastic_MITdiag_postResults + +contains + + +!-------------------------------------------------------------------------------------------------- +!> @brief module initialization +!> @details reads in material parameters, allocates arrays, and does sanity checks +!-------------------------------------------------------------------------------------------------- +subroutine plastic_MITdiag_init(fileUnit) + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) + use debug, only: & + debug_level, & + debug_constitutive, & + debug_levelBasic + use numerics, only: & + numerics_integrator + use math, only: & + math_Mandel3333to66, & + math_Voigt66to3333 + use IO, only: & + IO_read, & + IO_lc, & + IO_getTag, & + IO_isBlank, & + IO_stringPos, & + IO_stringValue, & + IO_intValue, & + IO_floatValue, & + IO_error, & + IO_timeStamp, & + IO_EOF, & + IO_warning + use lattice + use material, only: & + phase_plasticity, & + phase_plasticityInstance, & + phase_Noutput, & + PLASTICITY_MITdiag_label, & + PLASTICITY_MITdiag_ID, & + material_phase, & + plasticState, & + MATERIAL_partPhase + + use lattice + + implicit none + integer(pInt), intent(in) :: fileUnit + + + integer(pInt), allocatable, dimension(:) :: chunkPos + integer(pInt) :: & + f,i,j,o,k, & + phase, & + instance, & + maxNinstance, & + mySize = 0_pInt, & + sizeDotState, & + sizeState, & + sizeDeltaState, & + totalNslip, & + index_myFamily, & + index_otherFamily, & + Nchunks_SlipSlip = 0_pInt, & + Nchunks_SlipFamilies = 0_pInt + character(len=65536) :: & + tag = '', & + line = '', & + extmsg = '' + real(pReal), dimension(:), allocatable :: tempPerSlip + character(len=64) :: & + outputtag = '' + integer(pInt) :: NipcMyPhase + + write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_MITdiag_label//' init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() +#include "compilation_info.f90" + + maxNinstance = int(count(phase_plasticity == PLASTICITY_MITdiag_ID),pInt) + if (maxNinstance == 0_pInt) return + + if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & + write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance + + allocate(plastic_MITdiag_sizePostResults(maxNinstance), source=0_pInt) + allocate(plastic_MITdiag_sizePostResult(maxval(phase_Noutput), maxNinstance),source=0_pInt) + allocate(plastic_MITdiag_output(maxval(phase_Noutput), maxNinstance)) + plastic_MITdiag_output = '' + allocate(plastic_MITdiag_Noutput(maxNinstance), source=0_pInt) + allocate(plastic_MITdiag_totalNslip(maxNinstance), source=0_pInt) + allocate(param(maxNinstance)) ! one container of parameters per instance + + rewind(fileUnit) + phase = 0_pInt + do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partPhase) ! wind forward to + line = IO_read(fileUnit) + enddo + + parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part + line = IO_read(fileUnit) + if (IO_isBlank(line)) cycle ! skip empty lines + if (IO_getTag(line,'<','>') /= '') then ! stop at next part + line = IO_read(fileUnit, .true.) ! reset IO_read + exit + endif + if (IO_getTag(line,'[',']') /= '') then ! next section + phase = phase + 1_pInt ! advance section counter + if (phase_plasticity(phase) == PLASTICITY_MITdiag_ID) then + instance = phase_plasticityInstance(phase) ! count instances of my constitutive law + allocate(param(instance)%outputID(phase_Noutput(phase))) ! allocate space for IDs of every requested output + Nchunks_SlipFamilies = count(lattice_NslipSystem(:,phase) > 0_pInt) ! maximum number of slip families according to lattice type of current phase + Nchunks_SlipSlip = maxval(lattice_interactionSlipSlip(:,:,phase)) + if(allocated(tempPerSlip)) deallocate(tempPerSlip) + allocate(tempPerSlip(Nchunks_SlipFamilies)) + endif + cycle ! skip to next line + endif + + if (phase > 0_pInt) then; if (phase_plasticity(phase) == PLASTICITY_MITdiag_ID) then ! one of my phases. Do not short-circuit here (.and. between if-statements), it's not safe in Fortran + instance = phase_plasticityInstance(phase) ! which instance of my plasticity is present phase + chunkPos = IO_stringPos(line) + tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key + extmsg = trim(tag)//' ('//PLASTICITY_MITdiag_label//')' ! prepare error message identifier + + select case(tag) + case ('(output)') + outputtag = IO_lc(IO_stringValue(line,chunkPos,2_pInt)) + select case(outputtag) + case ('resistance') + plastic_MITdiag_Noutput(instance) = plastic_MITdiag_Noutput(instance) + 1_pInt + param(instance)%outputID(plastic_MITdiag_Noutput(instance)) = resistance_ID + plastic_MITdiag_output(plastic_MITdiag_Noutput(instance),instance) = outputtag + case ('accumulatedshear') + plastic_MITdiag_Noutput(instance) = plastic_MITdiag_Noutput(instance) + 1_pInt + param(instance)%outputID(plastic_MITdiag_Noutput(instance)) = accumulatedshear_ID + plastic_MITdiag_output(plastic_MITdiag_Noutput(instance),instance) = outputtag + case ('shearrate') + plastic_MITdiag_Noutput(instance) = plastic_MITdiag_Noutput(instance) + 1_pInt + param(instance)%outputID(plastic_MITdiag_Noutput(instance)) = shearrate_ID + plastic_MITdiag_output(plastic_MITdiag_Noutput(instance),instance) = outputtag + case ('resolvedstress') + plastic_MITdiag_Noutput(instance) = plastic_MITdiag_Noutput(instance) + 1_pInt + param(instance)%outputID(plastic_MITdiag_Noutput(instance)) = resolvedstress_ID + plastic_MITdiag_output(plastic_MITdiag_Noutput(instance),instance) = outputtag + + end select + +!-------------------------------------------------------------------------------------------------- +! parameters depending on number of slip families + case ('nslip') + if (chunkPos(1) < Nchunks_SlipFamilies + 1_pInt) & + call IO_warning(50_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_MITdiag_label//')') + if (chunkPos(1) > Nchunks_SlipFamilies + 1_pInt) & + call IO_error(150_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_MITdiag_label//')') + Nchunks_SlipFamilies = chunkPos(1) - 1_pInt ! user specified number of (possibly) active slip families (e.g. 6 0 6 --> 3) + do j = 1_pInt, Nchunks_SlipFamilies + param(instance)%n_slip(j) = IO_intValue(line,chunkPos,1_pInt+j) + enddo + case ('b','g0','rho0','rhosat','gammasat') + tempPerSlip = 0.0_pReal + do j = 1_pInt, Nchunks_SlipFamilies + if (param(instance)%n_slip(j) > 0_pInt) & + tempPerSlip(j) = IO_floatValue(line,chunkPos,1_pInt+j) + enddo + select case(tag) + case ('b') + param(instance)%b(1:Nchunks_SlipFamilies) = tempPerSlip(1:Nchunks_SlipFamilies) + case ('g0') + param(instance)%g0(1:Nchunks_SlipFamilies) = tempPerSlip(1:Nchunks_SlipFamilies) + case ('rho0') + param(instance)%rho0(1:Nchunks_SlipFamilies) = tempPerSlip(1:Nchunks_SlipFamilies) + case ('rhosat') + param(instance)%rhosat(1:Nchunks_SlipFamilies) = tempPerSlip(1:Nchunks_SlipFamilies) + case ('gammasat') + param(instance)%gammasat(1:Nchunks_SlipFamilies) = tempPerSlip(1:Nchunks_SlipFamilies) + end select + + case ('interaction_slipslip') + if (chunkPos(1) < 1_pInt + Nchunks_SlipSlip) & + call IO_warning(52_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_MITDIAG_label//')') + do j = 1_pInt, Nchunks_SlipSlip + param(instance)%interaction(j) = IO_floatValue(line,chunkPos,1_pInt+j) + enddo + + case ('gdot0') + param(instance)%gdot0 = IO_floatValue(line,chunkPos,2_pInt) + if (param(instance)%gdot0 <= 0.0_pReal) call IO_error(211_pInt,ext_msg=extmsg) + + case ('n') + param(instance)%n = IO_floatValue(line,chunkPos,2_pInt) + if (param(instance)%n <= 0.0_pReal) call IO_error(211_pInt,ext_msg=extmsg) + + case ('atol_resistance') + param(instance)%aTolResistance = IO_floatValue(line,chunkPos,2_pInt) + if (param(instance)%aTolResistance <= 0.0_pReal) call IO_error(211_pInt,ext_msg=extmsg) + + case ('atol_shear') + param(instance)%aTolShear = IO_floatValue(line,chunkPos,2_pInt) + + case default + + end select + endif; endif + enddo parsingFile + + allocate(state(maxNinstance)) ! internal state aliases + allocate(state0(maxNinstance)) + allocate(dotState(maxNinstance)) + allocate(stateAbsTol(maxNinstance)) + + initializeInstances: do phase = 1_pInt, size(phase_plasticity) ! loop over every plasticity + myPhase: if (phase_plasticity(phase) == PLASTICITY_MITdiag_ID) then ! isolate instances of own constitutive description + NipcMyPhase = count(material_phase == phase) ! number of own material points (including point components ipc) + instance = phase_plasticityInstance(phase) +!-------------------------------------------------------------------------------------------------- +! sanity checks + if (param(instance)%aTolShear <= 0.0_pReal) & + param(instance)%aTolShear = 1.0e-6_pReal ! default absolute tolerance 1e-6 + param(instance)%n_slip(1:lattice_maxNslipFamily) = & + min(lattice_NslipSystem(1:lattice_maxNslipFamily,phase),& ! limit active slip systems per family to min of available and requested + param(instance)%n_slip(1:lattice_maxNslipFamily)) + plastic_MITdiag_totalNslip(instance) = sum(param(instance)%n_slip(1:lattice_maxNslipFamily)) ! how many slip systems altogether + totalNslip = plastic_MITdiag_totalNslip(instance) + + allocate(plastic_MITdiag_hardeningMatrix(maxval(plastic_MITdiag_totalNslip),& ! slip resistance from slip activity + maxval(plastic_MITdiag_totalNslip),& + maxNinstance), source=0.0_pReal) + +!-------------------------------------------------------------------------------------------------- +! Determine size of postResults array + outputsLoop: do o = 1_pInt,plastic_MITdiag_Noutput(instance) + select case(param(instance)%outputID(o)) + case(resistance_ID, & + accumulatedshear_ID, & + shearrate_ID, & + resolvedstress_ID & + ) + mySize = plastic_MITdiag_totalNslip(instance) + case default + end select + + outputFound: if (mySize > 0_pInt) then + plastic_MITdiag_sizePostResult(o,instance) = mySize + plastic_MITdiag_sizePostResults(instance) = & + plastic_MITdiag_sizePostResults(instance) + mySize + endif outputFound + enddo outputsLoop + +!-------------------------------------------------------------------------------------------------- +! allocate state arrays + sizeState = totalNslip & ! g (resistance) + + totalNslip ! accshear + sizeDotState = sizeState ! both evolve + sizeDeltaState = 0_pInt ! no sudden jumps in state + plasticState(phase)%sizeState = sizeState + plasticState(phase)%sizeDotState = sizeDotState + plasticState(phase)%sizeDeltaState = sizeDeltaState + plasticState(phase)%sizePostResults = plastic_MITdiag_sizePostResults(instance) + plasticState(phase)%nSlip = totalNslip + plasticState(phase)%nTwin = 0 + plasticState(phase)%nTrans= 0 + allocate(plasticState(phase)%aTolState ( sizeState)) + + allocate(plasticState(phase)%state0 ( sizeState,NipcMyPhase),source=0.0_pReal) + + allocate(plasticState(phase)%partionedState0 ( sizeState,NipcMyPhase),source=0.0_pReal) + allocate(plasticState(phase)%subState0 ( sizeState,NipcMyPhase),source=0.0_pReal) + allocate(plasticState(phase)%state ( sizeState,NipcMyPhase),source=0.0_pReal) + allocate(plasticState(phase)%dotState (sizeDotState,NipcMyPhase),source=0.0_pReal) + allocate(plasticState(phase)%deltaState (sizeDeltaState,NipcMyPhase),source=0.0_pReal) + if (any(numerics_integrator == 1_pInt)) then + allocate(plasticState(phase)%previousDotState (sizeDotState,NipcMyPhase),source=0.0_pReal) + allocate(plasticState(phase)%previousDotState2(sizeDotState,NipcMyPhase),source=0.0_pReal) + endif + if (any(numerics_integrator == 4_pInt)) & + allocate(plasticState(phase)%RK4dotState (sizeDotState,NipcMyPhase),source=0.0_pReal) + if (any(numerics_integrator == 5_pInt)) & + allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NipcMyPhase),source=0.0_pReal) + + do f = 1_pInt,lattice_maxNslipFamily ! >>> interaction slip -- X + index_myFamily = sum(param(instance)%n_slip(1:f-1_pInt)) + do j = 1_pInt,param(instance)%n_slip(f) ! loop over (active) systems in my family (slip) + do o = 1_pInt,lattice_maxNslipFamily + index_otherFamily = sum(param(instance)%n_slip(1:o-1_pInt)) + do k = 1_pInt,param(instance)%n_slip(o) ! loop over (active) systems in other family (slip) + plastic_MITdiag_hardeningMatrix(index_myFamily+j, index_otherFamily+k, instance) = & + param(instance)%interaction(lattice_interactionSlipSlip( & + sum(lattice_NslipSystem(1:f-1,phase))+j, & + sum(lattice_NslipSystem(1:o-1,phase))+k, & + phase)) + enddo + enddo + enddo + enddo + +!-------------------------------------------------------------------------------------------------- +! globally required state aliases + plasticState(phase)%slipRate => plasticState(phase)%dotState(totalNslip+1:2*totalNslip,1:NipcMyPhase) + plasticState(phase)%accumulatedSlip => plasticState(phase)%state (totalNslip+1:2*totalNslip,1:NipcMyPhase) + +!-------------------------------------------------------------------------------------------------- +! locally defined state aliases + state(instance)%resistance => plasticState(phase)%state (1:totalNslip,1:NipcMyPhase) + state0(instance)%resistance => plasticState(phase)%state0 (1:totalNslip,1:NipcMyPhase) + dotState(instance)%resistance => plasticState(phase)%dotState (1:totalNslip,1:NipcMyPhase) + stateAbsTol(instance)%resistance => plasticState(phase)%aTolState(1) + + state(instance)%accumulatedShear => plasticState(phase)%state (totalNslip+1:2*totalNslip,1:NipcMyPhase) + state0(instance)%accumulatedShear => plasticState(phase)%state0 (totalNslip+1:2*totalNslip,1:NipcMyPhase) + dotState(instance)%accumulatedShear => plasticState(phase)%dotState (totalNslip+1:2*totalNslip,1:NipcMyPhase) + stateAbsTol(instance)%accumulatedShear => plasticState(phase)%aTolState(2) + +!-------------------------------------------------------------------------------------------------- +! init state + state0(instance)%resistance = param(instance)%g0(1) + state0(instance)%accumulatedShear = 1.0e-150_pReal + +!-------------------------------------------------------------------------------------------------- +! init absolute state tolerances + stateAbsTol(instance)%resistance = param(instance)%aTolResistance + stateAbsTol(instance)%accumulatedShear = param(instance)%aTolShear + + endif myPhase + enddo initializeInstances + +end subroutine plastic_MITdiag_init + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates plastic velocity gradient and its tangent +!-------------------------------------------------------------------------------------------------- +subroutine plastic_MITdiag_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el) + use prec, only: & + dNeq0 + use debug, only: & + debug_level, & + debug_constitutive, & + debug_levelBasic, & + debug_levelExtensive, & + debug_levelSelective, & + debug_e, & + debug_i, & + debug_g + use math, only: & + math_mul6x6, & + math_Mandel6to33, & + math_Plain3333to99, & + math_deviatoric33, & + math_mul33xx33, & + math_transpose33 + use lattice, only: & + lattice_Sslip, & + lattice_Sslip_v, & + lattice_maxNslipFamily, & + lattice_NslipSystem + use material, only: & + phaseAt, phasememberAt, & + material_phase, & + phase_plasticityInstance + + implicit none + real(pReal), dimension(3,3), intent(out) :: & + Lp !< plastic velocity gradient + real(pReal), dimension(9,9), intent(out) :: & + dLp_dTstar99 !< derivative of Lp with respect to 2nd Piola Kirchhoff stress + + real(pReal), dimension(6), intent(in) :: & + Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation + + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + + integer(pInt) :: & + instance, & + index_myFamily, & + f,i,j,k,l,m,n, & + of, & + ph + + real(pReal) :: & + tau_slip, & + gdot_slip, & + dgdot_dtauslip + + real(pReal), dimension(3,3,3,3) :: & + dLp_dTstar3333 !< derivative of Lp with respect to Tstar as 4th order tensor + + of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember + ph = phaseAt(ipc,ip,el) + instance = phase_plasticityInstance(material_phase(ipc,ip,el)) + + gdot_slip = 0.0_pReal + Lp = 0.0_pReal + dLp_dTstar3333 = 0.0_pReal + dLp_dTstar99 = 0.0_pReal + +!-------------------------------------------------------------------------------------------------- +! Slip part + j = 0_pInt + slipFamilies: do f = 1_pInt,lattice_maxNslipFamily + index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family + slipSystems: do i = 1_pInt,param(instance)%n_slip(f) + j = j+1_pInt + + ! Calculation of Lp + tau_slip = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph)) + + if (dNeq0(abs(tau_slip))) then + gdot_slip = param(instance)%gdot0 * & + ((abs(tau_slip)/(state(instance)%resistance(j,of))) & + **param(instance)%n)*sign(1.0_pReal,tau_slip) + + Lp = Lp + gdot_slip*lattice_Sslip(1:3,1:3,1,index_myFamily+i,ph) + end if + + ! Calculation of the tangent of Lp + if (dNeq0(abs(tau_slip))) then + dgdot_dtauslip = gdot_slip * param(instance)%n / tau_slip + forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & + dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + & + dgdot_dtauslip * lattice_Sslip(k,l,1,index_myFamily+i,ph)* & + lattice_Sslip(m,n,1,index_myFamily+i,ph) + endif + + enddo slipSystems + enddo slipFamilies + + dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333) + + if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0_pInt & + .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & + .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then + write(6,'(/,a,/,3(12x,3(f12.4,1x)/))') '<< CONST isotropic >> Lp', & + math_transpose33(Lp(1:3,1:3)) + write(6,'(/,a,/,f12.5)') '<< CONST phenopowerlaw >> gdot', gdot_slip + end if + +end subroutine plastic_MITdiag_LpAndItsTangent + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates the rate of change of microstructure +!-------------------------------------------------------------------------------------------------- +subroutine plastic_MITdiag_dotState(Tstar_v,ipc,ip,el) + use prec, only: & + dNeq0 + use math, only: & + pi + use material, only: & + phaseAt, phasememberAt, & + material_phase, & + phase_plasticityInstance + use lattice, only: & + lattice_mu, & + lattice_Sslip_v, & + lattice_maxNslipFamily, & + lattice_NslipSystem + + implicit none + real(pReal), dimension(6), intent(in):: & + Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + + integer(pInt) :: & + instance,ph, & + nSlip, & + f,i,j, & + index_myFamily, & + of + real(pReal) :: & + tau_slip, & + forresthardening + + real(pReal), dimension(plastic_MITdiag_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & + gdot_slip, & + tau_c, & + gamma_c, & + h, rho + + of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember + ph = phaseAt(ipc,ip,el) + instance = phase_plasticityInstance(material_phase(ipc,ip,el)) + nSlip = plastic_MITdiag_totalNslip(instance) + +!-------------------------------------------------------------------------------------------------- +! calculate left and right vectors and calculate dot gammas + gdot_slip = 0.0_pReal + tau_slip = 0.0_pReal + rho = 0.0_pReal + + j = 0_pInt + slipFamilies1: do f = 1_pInt,lattice_maxNslipFamily + index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family + slipSystems1: do i = 1_pInt,param(instance)%n_slip(f) + j = j+1_pInt + + tau_slip = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph)) + + if (dNeq0(abs(tau_slip))) then + gdot_slip(j) = param(instance)%gdot0 * & + ((abs(tau_slip)/(state(instance)%resistance(j,of))) & + **param(instance)%n * sign(1.0_pReal,tau_slip)) + + + endif + rho(j) = param(instance)%rhosat(f)* & + (1.0_pReal - (1.0_pReal - (param(instance)%rho0(f) / param(instance)%rhosat(f))) * & + exp(-state(instance)%accumulatedShear(j,of) / param(instance)%gammasat(f))) + + dotState(ph)%accumulatedShear(j,of) = abs(gdot_slip(j)) + + enddo slipSystems1 + enddo slipFamilies1 + + tau_c = 0.0_pReal + gamma_c = 0.0_pReal + h = 0.0_pReal + + + j = 0_pInt + slipFamilies2: do f = 1_pInt,lattice_maxNslipFamily + index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family + slipSystems2: do i = 1_pInt,param(instance)%n_slip(f) + j = j+1_pInt + + forresthardening = & + dot_product(plastic_MITdiag_hardeningMatrix(j,1:nSlip,instance), rho) + tau_c(j) = 0.3_pReal * lattice_mu(ph) * param(instance)%b(f) * & + sqrt(pi * forresthardening) + gamma_c(j) = param(instance)%b(f) * rho(j) / & + (2.0_pReal * sqrt(forresthardening)) + h(j) = (tau_c(j) / gamma_c(j)) * (state(instance)%resistance(j,of) / tau_c(j))**3.0_pReal * & + (cosh((tau_c(j) / state(instance)%resistance(j,of))**2.0_pReal) - 1.0_pReal) + dotState(ph)%resistance(j,of) = h(j) * abs(gdot_slip(j)) + + enddo slipSystems2 + enddo slipFamilies2 + +end subroutine plastic_MITdiag_dotState + +!-------------------------------------------------------------------------------------------------- +!> @brief return array of constitutive results +!-------------------------------------------------------------------------------------------------- +function plastic_MITdiag_postResults(Tstar_v,ipc,ip,el) + use math, only: & + math_mul6x6 + use material, only: & + material_phase, & + phaseAt, phasememberAt, & + phase_plasticityInstance + use lattice, only: & + lattice_NslipSystem, & + lattice_maxNslipFamily, & + lattice_Sslip_v + + implicit none + real(pReal), dimension(6), intent(in) :: & + Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), dimension(plastic_MITdiag_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & + plastic_MITdiag_postResults + + integer(pInt) :: & + instance,ph, of, & + nSlip, & + o,f,i,c,j, & + index_myFamily + real(pReal) :: & + tau_slip + + of = phasememberAt(ipc,ip,el) + ph = phaseAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember + instance = phase_plasticityInstance(material_phase(ipc,ip,el)) + + nSlip = plastic_MITdiag_totalNslip(instance) + + c = 0_pInt + plastic_MITdiag_postResults = 0.0_pReal + + outputsLoop: do o = 1_pInt,plastic_MITdiag_Noutput(instance) + select case(param(instance)%outputID(o)) + case (resistance_ID) + plastic_MITdiag_postResults(c+1_pInt:c+nSlip) = state(instance)%resistance(:,of) + c = c + nSlip + + case (accumulatedshear_ID) + plastic_MITdiag_postResults(c+1_pInt:c+nSlip) = state(instance)%accumulatedShear(:,of) + c = c + nSlip + + case (shearrate_ID) + j = 0_pInt + slipFamilies1: do f = 1_pInt,lattice_maxNslipFamily + index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family + slipSystems1: do i = 1_pInt,param(instance)%n_slip(f) + j = j + 1_pInt + tau_slip = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph)) + plastic_MITdiag_postResults(c+j) = param(instance)%gdot0 * & + ((abs(tau_slip) / state(instance)%resistance(j,of)) ** param(instance)%n & + *sign(1.0_pReal,tau_slip)) + enddo slipSystems1 + enddo slipFamilies1 + c = c + nSlip + + case (resolvedstress_ID) + j = 0_pInt + slipFamilies2: do f = 1_pInt,lattice_maxNslipFamily + index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family + slipSystems2: do i = 1_pInt,param(instance)%n_slip(f) + j = j + 1_pInt + plastic_MITdiag_postResults(c+j) = & + dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph)) + enddo slipSystems2 + enddo slipFamilies2 + c = c + nSlip + + end select + enddo outputsLoop + +end function plastic_MITdiag_postResults + + +end module plastic_MITdiag