From dd68374afd452c2fe3da1a8c2970b87adf93de85 Mon Sep 17 00:00:00 2001 From: Tias Maiti Date: Sun, 11 Jun 2017 18:47:21 -0400 Subject: [PATCH] moved new constitutive law to new branch for further testing --- src/CMakeLists.txt | 1 - src/constitutive.f90 | 27 +- src/material.f90 | 26 +- src/plastic_MITdiag.f90 | 716 ---------------------------------------- 4 files changed, 15 insertions(+), 755 deletions(-) delete mode 100644 src/plastic_MITdiag.f90 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 6b32e39ec..eb1448028 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -74,7 +74,6 @@ 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 99594a846..de8f61c2a 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -70,7 +70,6 @@ subroutine constitutive_init() PLASTICITY_none_ID, & PLASTICITY_isotropic_ID, & PLASTICITY_phenopowerlaw_ID, & - PLASTICITY_mitdiag_ID ,& PLASTICITY_phenoplus_ID, & PLASTICITY_dislotwin_ID, & PLASTICITY_disloucla_ID, & @@ -94,7 +93,6 @@ subroutine constitutive_init() PLASTICITY_NONE_label, & PLASTICITY_ISOTROPIC_label, & PLASTICITY_PHENOPOWERLAW_label, & - PLASTICITY_MITDIAG_label, & PLASTICITY_PHENOPLUS_label, & PLASTICITY_DISLOTWIN_label, & PLASTICITY_DISLOUCLA_label, & @@ -115,7 +113,6 @@ subroutine constitutive_init() use plastic_none use plastic_isotropic use plastic_phenopowerlaw - use plastic_mitdiag use plastic_phenoplus use plastic_dislotwin use plastic_disloucla @@ -161,7 +158,6 @@ 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) @@ -211,7 +207,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 @@ -222,11 +218,6 @@ 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 @@ -510,7 +501,6 @@ 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, & @@ -520,8 +510,6 @@ 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: & @@ -572,8 +560,6 @@ 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 @@ -898,7 +884,6 @@ 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, & @@ -912,8 +897,6 @@ 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: & @@ -967,8 +950,6 @@ 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 @@ -1112,7 +1093,6 @@ 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, & @@ -1128,8 +1108,6 @@ 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: & @@ -1182,9 +1160,6 @@ 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 1acff4f93..a77c4871a 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -25,7 +25,6 @@ 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', & @@ -74,7 +73,6 @@ module material enumerator :: PLASTICITY_undefined_ID, & PLASTICITY_none_ID, & PLASTICITY_isotropic_ID, & - PLASTICITY_mitdiag_ID, & PLASTICITY_phenopowerlaw_ID, & PLASTICITY_phenoplus_ID, & PLASTICITY_dislotwin_ID, & @@ -314,7 +312,6 @@ module material PLASTICITY_none_ID, & PLASTICITY_isotropic_ID, & PLASTICITY_phenopowerlaw_ID, & - PLASTICITY_mitdiag_ID, & PLASTICITY_phenoplus_ID, & PLASTICITY_dislotwin_ID, & PLASTICITY_disloucla_ID, & @@ -674,7 +671,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 @@ -732,6 +729,8 @@ 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, & @@ -741,7 +740,6 @@ 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) :: & @@ -761,7 +759,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='Microstructure in geometry > Sections in material.config') + call IO_error(155_pInt,ext_msg='More microstructures in geometry than 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 @@ -804,7 +802,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') @@ -819,6 +817,12 @@ 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 @@ -975,19 +979,17 @@ 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 deleted file mode 100644 index 28d9bc3d7..000000000 --- a/src/plastic_MITdiag.f90 +++ /dev/null @@ -1,716 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @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