diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 713aab5d7..8a12bef76 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -19,7 +19,7 @@ module CPFEM use HDF5_utilities use results use lattice - use constitutive + use phase implicit none private diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index e696858cf..2b32a0cbb 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -19,7 +19,7 @@ module CPFEM2 use discretization use HDF5_utilities use homogenization - use constitutive + use phase #if defined(Mesh) use FEM_quadrature use discretization_mesh diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index a20fb1cee..6480f7382 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -7,7 +7,6 @@ #include "IO.f90" #include "YAML_types.f90" #include "YAML_parse.f90" -#include "future.f90" #include "config.f90" #include "LAPACK_interface.f90" #include "math.f90" @@ -17,38 +16,38 @@ #include "results.f90" #include "geometry_plastic_nonlocal.f90" #include "discretization.f90" -#ifdef Marc4DAMASK #include "marc/discretization_marc.f90" -#endif #include "material.f90" #include "lattice.f90" -#include "constitutive.f90" -#include "constitutive_mech.f90" -#include "constitutive_plastic_none.f90" -#include "constitutive_plastic_isotropic.f90" -#include "constitutive_plastic_phenopowerlaw.f90" -#include "constitutive_plastic_kinehardening.f90" -#include "constitutive_plastic_dislotwin.f90" -#include "constitutive_plastic_disloTungsten.f90" -#include "constitutive_plastic_nonlocal.f90" -#include "constitutive_thermal.f90" -#include "constitutive_thermal_dissipation.f90" -#include "constitutive_thermal_externalheat.f90" -#include "kinematics_thermal_expansion.f90" -#include "constitutive_damage.f90" -#include "source_damage_isoBrittle.f90" -#include "source_damage_isoDuctile.f90" -#include "source_damage_anisoBrittle.f90" -#include "source_damage_anisoDuctile.f90" -#include "kinematics_cleavage_opening.f90" -#include "kinematics_slipplane_opening.f90" +#include "phase.f90" +#include "phase_mechanics.f90" +#include "phase_mechanics_plastic.f90" +#include "phase_mechanics_plastic_none.f90" +#include "phase_mechanics_plastic_isotropic.f90" +#include "phase_mechanics_plastic_phenopowerlaw.f90" +#include "phase_mechanics_plastic_kinehardening.f90" +#include "phase_mechanics_plastic_dislotwin.f90" +#include "phase_mechanics_plastic_dislotungsten.f90" +#include "phase_mechanics_plastic_nonlocal.f90" +#include "phase_mechanics_eigendeformation.f90" +#include "phase_mechanics_eigendeformation_cleavageopening.f90" +#include "phase_mechanics_eigendeformation_slipplaneopening.f90" +#include "phase_mechanics_eigendeformation_thermalexpansion.f90" +#include "phase_thermal.f90" +#include "phase_thermal_dissipation.f90" +#include "phase_thermal_externalheat.f90" +#include "phase_damage.f90" +#include "phase_damage_isobrittle.f90" +#include "phase_damage_isoductile.f90" +#include "phase_damage_anisobrittle.f90" +#include "phase_damage_anisoductile.f90" #include "damage_none.f90" #include "damage_nonlocal.f90" #include "homogenization.f90" -#include "homogenization_mech.f90" -#include "homogenization_mech_none.f90" -#include "homogenization_mech_isostrain.f90" -#include "homogenization_mech_RGC.f90" +#include "homogenization_mechanics.f90" +#include "homogenization_mechanics_none.f90" +#include "homogenization_mechanics_isostrain.f90" +#include "homogenization_mechanics_RGC.f90" #include "homogenization_thermal.f90" #include "homogenization_damage.f90" #include "CPFEM.f90" diff --git a/src/constitutive_thermal_externalheat.f90 b/src/constitutive_thermal_externalheat.f90 deleted file mode 100644 index 16b2bdb65..000000000 --- a/src/constitutive_thermal_externalheat.f90 +++ /dev/null @@ -1,137 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH -!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH -!> @author Philip Eisenlohr, Michigan State University -!> @brief material subroutine for variable heat source -!-------------------------------------------------------------------------------------------------- -submodule(constitutive:constitutive_thermal) source_externalheat - - - integer, dimension(:), allocatable :: & - source_thermal_externalheat_offset, & !< which source is my current thermal dissipation mechanism? - source_thermal_externalheat_instance !< instance of thermal dissipation source mechanism - - type :: tParameters !< container type for internal constitutive parameters - real(pReal), dimension(:), allocatable :: & - t_n, & - f_T - integer :: & - nIntervals - end type tParameters - - type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances) - - -contains - - -!-------------------------------------------------------------------------------------------------- -!> @brief module initialization -!> @details reads in material parameters, allocates arrays, and does sanity checks -!-------------------------------------------------------------------------------------------------- -module function source_thermal_externalheat_init(source_length) result(mySources) - - integer, intent(in) :: source_length - logical, dimension(:,:), allocatable :: mySources - - class(tNode), pointer :: & - phases, & - phase, & - sources, thermal, & - src - integer :: Ninstances,sourceOffset,Nconstituents,p - - print'(/,a)', ' <<<+- thermal_externalheat init -+>>>' - - mySources = thermal_active('externalheat',source_length) - - Ninstances = count(mySources) - print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) - if(Ninstances == 0) return - - phases => config_material%get('phase') - allocate(param(Ninstances)) - allocate(source_thermal_externalheat_offset (phases%length), source=0) - allocate(source_thermal_externalheat_instance(phases%length), source=0) - - do p = 1, phases%length - phase => phases%get(p) - if(any(mySources(:,p))) source_thermal_externalheat_instance(p) = count(mySources(:,1:p)) - if(count(mySources(:,p)) == 0) cycle - thermal => phase%get('thermal') - sources => thermal%get('source') - do sourceOffset = 1, sources%length - if(mySources(sourceOffset,p)) then - source_thermal_externalheat_offset(p) = sourceOffset - associate(prm => param(source_thermal_externalheat_instance(p))) - src => sources%get(sourceOffset) - - prm%t_n = src%get_asFloats('t_n') - prm%nIntervals = size(prm%t_n) - 1 - - prm%f_T = src%get_asFloats('f_T',requiredSize = size(prm%t_n)) - - Nconstituents = count(material_phaseAt==p) * discretization_nIPs - call constitutive_allocateState(thermalState(p)%p(sourceOffset),Nconstituents,1,1,0) - end associate - endif - enddo - enddo - -end function source_thermal_externalheat_init - - -!-------------------------------------------------------------------------------------------------- -!> @brief rate of change of state -!> @details state only contains current time to linearly interpolate given heat powers -!-------------------------------------------------------------------------------------------------- -module subroutine source_thermal_externalheat_dotState(ph, me) - - integer, intent(in) :: & - ph, & - me - - integer :: & - sourceOffset - - sourceOffset = source_thermal_externalheat_offset(ph) - - thermalState(ph)%p(sourceOffset)%dotState(1,me) = 1.0_pReal ! state is current time - -end subroutine source_thermal_externalheat_dotState - - -!-------------------------------------------------------------------------------------------------- -!> @brief returns local heat generation rate -!-------------------------------------------------------------------------------------------------- -module subroutine thermal_externalheat_getRate(TDot, ph, me) - - integer, intent(in) :: & - ph, & - me - real(pReal), intent(out) :: & - TDot - - integer :: & - sourceOffset, interval - real(pReal) :: & - frac_time - - sourceOffset = source_thermal_externalheat_offset(ph) - - associate(prm => param(source_thermal_externalheat_instance(ph))) - do interval = 1, prm%nIntervals ! scan through all rate segments - frac_time = (thermalState(ph)%p(sourceOffset)%state(1,me) - prm%t_n(interval)) & - / (prm%t_n(interval+1) - prm%t_n(interval)) ! fractional time within segment - if ( (frac_time < 0.0_pReal .and. interval == 1) & - .or. (frac_time >= 1.0_pReal .and. interval == prm%nIntervals) & - .or. (frac_time >= 0.0_pReal .and. frac_time < 1.0_pReal) ) & - TDot = prm%f_T(interval ) * (1.0_pReal - frac_time) + & - prm%f_T(interval+1) * frac_time ! interpolate heat rate between segment boundaries... - ! ...or extrapolate if outside me bounds - enddo - end associate - -end subroutine thermal_externalheat_getRate - -end submodule source_externalheat diff --git a/src/damage_nonlocal.f90 b/src/damage_nonlocal.f90 index f566ebbeb..807231889 100644 --- a/src/damage_nonlocal.f90 +++ b/src/damage_nonlocal.f90 @@ -8,7 +8,7 @@ module damage_nonlocal use config use YAML_types use lattice - use constitutive + use phase use results implicit none diff --git a/src/future.f90 b/src/future.f90 deleted file mode 100644 index b7eb3fec9..000000000 --- a/src/future.f90 +++ /dev/null @@ -1,35 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH -!> @brief New fortran functions for compiler versions that do not support them -!-------------------------------------------------------------------------------------------------- -module future - use prec - - implicit none - public - -contains - -#if defined(__GFORTRAN__) && __GNUC__<9 || defined(__INTEL_COMPILER) && INTEL_COMPILER<1800 -!-------------------------------------------------------------------------------------------------- -!> @brief substitute for the findloc intrinsic (only for integer, dimension(:) at the moment) -!-------------------------------------------------------------------------------------------------- -function findloc(a,v) - - integer, intent(in), dimension(:) :: a - integer, intent(in) :: v - integer :: i,j - integer, allocatable, dimension(:) :: findloc - - allocate(findloc(count(a==v))) - j = 1 - do i = 1, size(a) - if (a(i)==v) then - findloc(j) = i - j = j + 1 - endif - enddo -end function findloc -#endif - -end module future diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 9b8c33c2e..2255e83b0 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -10,7 +10,7 @@ module homogenization use config use math use material - use constitutive + use phase use discretization use damage_none use damage_nonlocal @@ -145,8 +145,6 @@ module homogenization integer, intent(in) :: & ip, & !< integration point number el !< element number - integer :: & - co real(pReal) :: M end function damage_nonlocal_getMobility diff --git a/src/homogenization_damage.f90 b/src/homogenization_damage.f90 index 3115e6a6f..ba021a0e0 100644 --- a/src/homogenization_damage.f90 +++ b/src/homogenization_damage.f90 @@ -34,8 +34,8 @@ module subroutine damage_init() integer :: ho - print'(/,a)', ' <<<+- homogenization_damage init -+>>>' - + print'(/,a)', ' <<<+- homogenization:damage init -+>>>' + print'(/,a)', ' <<<+- homogenization:damage:isodamage init -+>>>' configHomogenizations => config_material%get('homogenization') allocate(param(configHomogenizations%length)) diff --git a/src/homogenization_mech.f90 b/src/homogenization_mechanics.f90 similarity index 98% rename from src/homogenization_mech.f90 rename to src/homogenization_mechanics.f90 index dd3e47c59..4eae859ee 100644 --- a/src/homogenization_mech.f90 +++ b/src/homogenization_mechanics.f90 @@ -2,7 +2,7 @@ !> @author Martin Diehl, KU Leuven !> @brief Partition F and homogenize P/dPdF !-------------------------------------------------------------------------------------------------- -submodule(homogenization) homogenization_mech +submodule(homogenization) mechanics interface @@ -86,7 +86,7 @@ module subroutine mech_init(num_homog) class(tNode), pointer :: & num_homogMech - print'(/,a)', ' <<<+- homogenization_mech init -+>>>' + print'(/,a)', ' <<<+- homogenization:mechanics init -+>>>' allocate(homogenization_dPdF(3,3,3,3,discretization_nIPs*discretization_Nelems), source=0.0_pReal) homogenization_F0 = spread(math_I3,3,discretization_nIPs*discretization_Nelems) ! initialize to identity @@ -253,4 +253,4 @@ module subroutine mech_results(group_base,h) end subroutine mech_results -end submodule homogenization_mech +end submodule mechanics diff --git a/src/homogenization_mech_RGC.f90 b/src/homogenization_mechanics_RGC.f90 similarity index 99% rename from src/homogenization_mech_RGC.f90 rename to src/homogenization_mechanics_RGC.f90 index 2c7a5a8cb..89a84314b 100644 --- a/src/homogenization_mech_RGC.f90 +++ b/src/homogenization_mechanics_RGC.f90 @@ -6,7 +6,7 @@ !> @brief Relaxed grain cluster (RGC) homogenization scheme !> N_constituents is defined as p x q x r (cluster) !-------------------------------------------------------------------------------------------------- -submodule(homogenization:homogenization_mech) homogenization_mech_RGC +submodule(homogenization:mechanics) RGC use rotations use lattice @@ -88,7 +88,7 @@ module subroutine mech_RGC_init(num_homogMech) homog, & homogMech - print'(/,a)', ' <<<+- homogenization_mech_rgc init -+>>>' + print'(/,a)', ' <<<+- homogenization:mechanics:RGC init -+>>>' Ninstances = count(homogenization_type == HOMOGENIZATION_RGC_ID) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) @@ -947,4 +947,4 @@ pure function interface1to4(iFace1D, nGDim) end function interface1to4 -end submodule homogenization_mech_RGC +end submodule RGC diff --git a/src/homogenization_mech_isostrain.f90 b/src/homogenization_mechanics_isostrain.f90 similarity index 96% rename from src/homogenization_mech_isostrain.f90 rename to src/homogenization_mechanics_isostrain.f90 index df8f5fc9d..9634b2a38 100644 --- a/src/homogenization_mech_isostrain.f90 +++ b/src/homogenization_mechanics_isostrain.f90 @@ -4,7 +4,7 @@ !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @brief Isostrain (full constraint Taylor assuption) homogenization scheme !-------------------------------------------------------------------------------------------------- -submodule(homogenization:homogenization_mech) homogenization_mech_isostrain +submodule(homogenization:mechanics) isostrain enum, bind(c); enumerator :: & parallel_ID, & @@ -37,7 +37,7 @@ module subroutine mech_isostrain_init homog, & homogMech - print'(/,a)', ' <<<+- homogenization_mech_isostrain init -+>>>' + print'(/,a)', ' <<<+- homogenization:mechanics:isostrain init -+>>>' Ninstances = count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) @@ -114,4 +114,4 @@ module subroutine mech_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dP end subroutine mech_isostrain_averageStressAndItsTangent -end submodule homogenization_mech_isostrain +end submodule isostrain diff --git a/src/homogenization_mech_none.f90 b/src/homogenization_mechanics_none.f90 similarity index 89% rename from src/homogenization_mech_none.f90 rename to src/homogenization_mechanics_none.f90 index 767dbf557..ebe2ea8f1 100644 --- a/src/homogenization_mech_none.f90 +++ b/src/homogenization_mechanics_none.f90 @@ -4,7 +4,7 @@ !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief dummy homogenization homogenization scheme for 1 constituent per material point !-------------------------------------------------------------------------------------------------- -submodule(homogenization:homogenization_mech) homogenization_mech_none +submodule(homogenization:mechanics) none contains @@ -18,7 +18,7 @@ module subroutine mech_none_init h, & Nmaterialpoints - print'(/,a)', ' <<<+- homogenization_mech_none init -+>>>' + print'(/,a)', ' <<<+- homogenization:mechanics:none init -+>>>' Ninstances = count(homogenization_type == HOMOGENIZATION_NONE_ID) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) @@ -38,4 +38,4 @@ module subroutine mech_none_init end subroutine mech_none_init -end submodule homogenization_mech_none +end submodule none diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90 index 6a12730ac..b7170e6ec 100644 --- a/src/homogenization_thermal.f90 +++ b/src/homogenization_thermal.f90 @@ -1,7 +1,7 @@ !-------------------------------------------------------------------------------------------------- !> @author Martin Diehl, KU Leuven !-------------------------------------------------------------------------------------------------- -submodule(homogenization) homogenization_thermal +submodule(homogenization) thermal use lattice @@ -34,7 +34,9 @@ module subroutine thermal_init() integer :: ho - print'(/,a)', ' <<<+- homogenization_thermal init -+>>>' + print'(/,a)', ' <<<+- homogenization:thermal init -+>>>' + print'(/,a)', ' <<<+- homogenization:thermal:isotemperature init -+>>>' + configHomogenizations => config_material%get('homogenization') @@ -225,15 +227,21 @@ module subroutine thermal_conduction_getSource(Tdot, ip,el) real(pReal), intent(out) :: & Tdot - integer :: & - homog + integer :: co, ho,ph,me + real(pReal) :: dot_T_temp - homog = material_homogenizationAt(el) - call constitutive_thermal_getRate(TDot, ip,el) + ho = material_homogenizationAt(el) + Tdot = 0.0_pReal + do co = 1, homogenization_Nconstituents(ho) + ph = material_phaseAt(co,el) + me = material_phasememberAt(co,ip,el) + call constitutive_thermal_getRate(dot_T_temp, ph,me) + Tdot = Tdot + dot_T_temp + enddo - Tdot = Tdot/real(homogenization_Nconstituents(homog),pReal) + Tdot = Tdot/real(homogenization_Nconstituents(ho),pReal) end subroutine thermal_conduction_getSource -end submodule homogenization_thermal +end submodule thermal diff --git a/src/constitutive.f90 b/src/phase.f90 similarity index 81% rename from src/constitutive.f90 rename to src/phase.f90 index 1412f5259..dd14e0831 100644 --- a/src/constitutive.f90 +++ b/src/phase.f90 @@ -3,7 +3,7 @@ !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @brief elasticity, plasticity, damage & thermal internal microstructure state !-------------------------------------------------------------------------------------------------- -module constitutive +module phase use prec use math use rotations @@ -15,18 +15,10 @@ module constitutive use discretization use parallelization use HDF5_utilities - use results implicit none private - enum, bind(c); enumerator :: & - KINEMATICS_UNDEFINED_ID ,& - KINEMATICS_CLEAVAGE_OPENING_ID, & - KINEMATICS_SLIPPLANE_OPENING_ID, & - KINEMATICS_THERMAL_EXPANSION_ID - end enum - type(rotation), dimension(:,:,:), allocatable :: & crystallite_orientation !< current orientation @@ -65,10 +57,6 @@ module constitutive type(tDebugOptions) :: debugCrystallite - - integer(kind(KINEMATICS_UNDEFINED_ID)), dimension(:,:), allocatable :: & - phase_kinematics !< active kinematic mechanisms of each phase - integer, dimension(:), allocatable, public :: & !< ToDo: should be protected (bug in Intel compiler) thermal_Nsources, & phase_Nsources, & !< number of source mechanisms active in each phase @@ -188,6 +176,11 @@ module constitutive real(pReal) :: T end function thermal_T + module function thermal_dot_T(ph,me) result(dot_T) + integer, intent(in) :: ph,me + real(pReal) :: dot_T + end function thermal_dot_T + module subroutine constitutive_mech_setF(F,co,ip,el) real(pReal), dimension(3,3), intent(in) :: F @@ -246,16 +239,12 @@ module constitutive dPhiDot_dPhi end subroutine constitutive_damage_getRateAndItsTangents - module subroutine constitutive_thermal_getRate(TDot, ip,el) - integer, intent(in) :: & - ip, & !< integration point number - el !< element number + module subroutine constitutive_thermal_getRate(TDot, ph,me) + integer, intent(in) :: ph, me real(pReal), intent(out) :: & TDot end subroutine constitutive_thermal_getRate - - module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e) integer, intent(in) :: & instance, & @@ -265,66 +254,31 @@ module constitutive orientation !< crystal orientation end subroutine plastic_nonlocal_updateCompatibility - - module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of) - real(pReal), dimension(3,3), intent(out) :: & - Li !< inleastic velocity gradient - real(pReal), dimension(3,3,3,3), intent(out) :: & - dLi_dMi !< derivative of Li with respect to Mandel stress - real(pReal), dimension(3,3), intent(in) :: & - Mi !< Mandel stress - integer, intent(in) :: & - instance, & - of - end subroutine plastic_isotropic_LiAndItsTangent - - module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el) - integer, intent(in) :: & - co, & !< grain number - ip, & !< integration point number - el !< element number - real(pReal), intent(in), dimension(3,3) :: & - S - real(pReal), intent(out), dimension(3,3) :: & - Ld !< damage velocity gradient - real(pReal), intent(out), dimension(3,3,3,3) :: & - dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) - end subroutine kinematics_cleavage_opening_LiAndItsTangent - - module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el) - integer, intent(in) :: & - co, & !< grain number - ip, & !< integration point number - el !< element number - real(pReal), intent(in), dimension(3,3) :: & - S - real(pReal), intent(out), dimension(3,3) :: & - Ld !< damage velocity gradient - real(pReal), intent(out), dimension(3,3,3,3) :: & - dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) - end subroutine kinematics_slipplane_opening_LiAndItsTangent - - module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ph,me) - integer, intent(in) :: ph, me - !< element number - real(pReal), intent(out), dimension(3,3) :: & - Li !< thermal velocity gradient - real(pReal), intent(out), dimension(3,3,3,3) :: & - dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero) - end subroutine kinematics_thermal_expansion_LiAndItsTangent - - module subroutine constitutive_plastic_dependentState(co,ip,el) + module subroutine plastic_dependentState(co,ip,el) integer, intent(in) :: & co, & !< component-ID of integration point ip, & !< integration point el !< element - end subroutine constitutive_plastic_dependentState + end subroutine plastic_dependentState end interface type(tDebugOptions) :: debugConstitutive +#if __INTEL_COMPILER >= 1900 + public :: & + prec, & + math, & + rotations, & + IO, & + config, & + material, & + results, & + lattice, & + discretization, & + HDF5_utilities +#endif public :: & constitutive_init, & @@ -336,7 +290,6 @@ module constitutive constitutive_forward, & constitutive_restore, & plastic_nonlocal_updateCompatibility, & - kinematics_active, & converged, & crystallite_init, & crystallite_stress, & @@ -353,15 +306,10 @@ module constitutive constitutive_mech_getP, & constitutive_mech_setF, & constitutive_mech_getF, & - constitutive_windForward, & - KINEMATICS_UNDEFINED_ID ,& - KINEMATICS_CLEAVAGE_OPENING_ID, & - KINEMATICS_SLIPPLANE_OPENING_ID, & - KINEMATICS_THERMAL_EXPANSION_ID + constitutive_windForward contains - !-------------------------------------------------------------------------------------------------- !> @brief Initialze constitutive models for individual physics !-------------------------------------------------------------------------------------------------- @@ -375,7 +323,7 @@ subroutine constitutive_init phases - print'(/,a)', ' <<<+- constitutive init -+>>>'; flush(IO_STDOUT) + print'(/,a)', ' <<<+- phase init -+>>>'; flush(IO_STDOUT) debug_constitutive => config_debug%get('constitutive', defaultVal=emptyList) debugConstitutive%basic = debug_constitutive%contains('basic') @@ -410,38 +358,6 @@ subroutine constitutive_init end subroutine constitutive_init - -!-------------------------------------------------------------------------------------------------- -!> @brief checks if a kinematic mechanism is active or not -!-------------------------------------------------------------------------------------------------- -function kinematics_active(kinematics_label,kinematics_length) result(active_kinematics) - - character(len=*), intent(in) :: kinematics_label !< name of kinematic mechanism - integer, intent(in) :: kinematics_length !< max. number of kinematics in system - logical, dimension(:,:), allocatable :: active_kinematics - - class(tNode), pointer :: & - phases, & - phase, & - kinematics, & - kinematics_type - integer :: p,k - - phases => config_material%get('phase') - allocate(active_kinematics(kinematics_length,phases%length), source = .false. ) - do p = 1, phases%length - phase => phases%get(p) - kinematics => phase%get('kinematics',defaultVal=emptyList) - do k = 1, kinematics%length - kinematics_type => kinematics%get(k) - if(kinematics_type%get_asString('type') == kinematics_label) active_kinematics(k,p) = .true. - enddo - enddo - - -end function kinematics_active - - !-------------------------------------------------------------------------------------------------- !> @brief Allocate the components of the state structure for a given phase !-------------------------------------------------------------------------------------------------- @@ -562,9 +478,7 @@ subroutine crystallite_init() class(tNode), pointer :: & num_crystallite, & debug_crystallite, & ! pointer to debug options for crystallite - phases, & - phase, & - mech + phases print'(/,a)', ' <<<+- crystallite init -+>>>' @@ -630,7 +544,7 @@ subroutine crystallite_init() ph = material_phaseAt(co,el) me = material_phaseMemberAt(co,ip,el) call crystallite_orientations(co,ip,el) - call constitutive_plastic_dependentState(co,ip,el) ! update dependent state variables to be consistent with basic states + call plastic_dependentState(co,ip,el) ! update dependent state variables to be consistent with basic states enddo enddo enddo @@ -788,4 +702,4 @@ subroutine constitutive_restartRead(fileHandle) end subroutine constitutive_restartRead -end module constitutive +end module phase diff --git a/src/constitutive_damage.f90 b/src/phase_damage.f90 similarity index 66% rename from src/constitutive_damage.f90 rename to src/phase_damage.f90 index 07bfbeeca..014347c1a 100644 --- a/src/constitutive_damage.f90 +++ b/src/phase_damage.f90 @@ -1,7 +1,7 @@ !---------------------------------------------------------------------------------------------------- !> @brief internal microstructure state for all damage sources and kinematics constitutive models !---------------------------------------------------------------------------------------------------- -submodule(constitutive) constitutive_damage +submodule(phase) damagee enum, bind(c); enumerator :: & DAMAGE_UNDEFINED_ID, & DAMAGE_ISOBRITTLE_ID, & @@ -22,35 +22,26 @@ submodule(constitutive) constitutive_damage interface - module function source_damage_anisoBrittle_init(source_length) result(mySources) - integer, intent(in) :: source_length - logical, dimension(:,:), allocatable :: mySources - end function source_damage_anisoBrittle_init + module function anisobrittle_init(source_length) result(mySources) + integer, intent(in) :: source_length + logical, dimension(:,:), allocatable :: mySources + end function anisobrittle_init - module function source_damage_anisoDuctile_init(source_length) result(mySources) - integer, intent(in) :: source_length - logical, dimension(:,:), allocatable :: mySources - end function source_damage_anisoDuctile_init + module function anisoductile_init(source_length) result(mySources) + integer, intent(in) :: source_length + logical, dimension(:,:), allocatable :: mySources + end function anisoductile_init - module function source_damage_isoBrittle_init(source_length) result(mySources) - integer, intent(in) :: source_length - logical, dimension(:,:), allocatable :: mySources - end function source_damage_isoBrittle_init + module function isobrittle_init(source_length) result(mySources) + integer, intent(in) :: source_length + logical, dimension(:,:), allocatable :: mySources + end function isobrittle_init - module function source_damage_isoDuctile_init(source_length) result(mySources) - integer, intent(in) :: source_length - logical, dimension(:,:), allocatable :: mySources - end function source_damage_isoDuctile_init + module function isoductile_init(source_length) result(mySources) + integer, intent(in) :: source_length + logical, dimension(:,:), allocatable :: mySources + end function isoductile_init - module function kinematics_cleavage_opening_init(kinematics_length) result(myKinematics) - integer, intent(in) :: kinematics_length - logical, dimension(:,:), allocatable :: myKinematics - end function kinematics_cleavage_opening_init - - module function kinematics_slipplane_opening_init(kinematics_length) result(myKinematics) - integer, intent(in) :: kinematics_length - logical, dimension(:,:), allocatable :: myKinematics - end function kinematics_slipplane_opening_init module subroutine source_damage_isoBrittle_deltaState(C, Fe, ph, me) integer, intent(in) :: ph,me @@ -61,93 +52,93 @@ submodule(constitutive) constitutive_damage end subroutine source_damage_isoBrittle_deltaState - module subroutine source_damage_anisoBrittle_dotState(S, co, ip, el) + module subroutine anisobrittle_dotState(S, co, ip, el) integer, intent(in) :: & co, & !< component-ID of integration point ip, & !< integration point el !< element real(pReal), intent(in), dimension(3,3) :: & S - end subroutine source_damage_anisoBrittle_dotState + end subroutine anisobrittle_dotState - module subroutine source_damage_anisoDuctile_dotState(co, ip, el) + module subroutine anisoductile_dotState(co, ip, el) integer, intent(in) :: & co, & !< component-ID of integration point ip, & !< integration point el !< element - end subroutine source_damage_anisoDuctile_dotState + end subroutine anisoductile_dotState - module subroutine source_damage_isoDuctile_dotState(co, ip, el) + module subroutine isoductile_dotState(co, ip, el) integer, intent(in) :: & co, & !< component-ID of integration point ip, & !< integration point el !< element - end subroutine source_damage_isoDuctile_dotState + end subroutine isoductile_dotState - module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - integer, intent(in) :: & - phase, & !< phase ID of element - constituent !< position of element within its phase instance - real(pReal), intent(in) :: & - phi !< damage parameter - real(pReal), intent(out) :: & - localphiDot, & - dLocalphiDot_dPhi - end subroutine source_damage_anisoBrittle_getRateAndItsTangent + module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) + integer, intent(in) :: & + phase, & !< phase ID of element + constituent !< position of element within its phase instance + real(pReal), intent(in) :: & + phi !< damage parameter + real(pReal), intent(out) :: & + localphiDot, & + dLocalphiDot_dPhi + end subroutine source_damage_anisoBrittle_getRateAndItsTangent - module subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - integer, intent(in) :: & - phase, & !< phase ID of element - constituent !< position of element within its phase instance - real(pReal), intent(in) :: & - phi !< damage parameter - real(pReal), intent(out) :: & - localphiDot, & - dLocalphiDot_dPhi - end subroutine source_damage_anisoDuctile_getRateAndItsTangent + module subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) + integer, intent(in) :: & + phase, & !< phase ID of element + constituent !< position of element within its phase instance + real(pReal), intent(in) :: & + phi !< damage parameter + real(pReal), intent(out) :: & + localphiDot, & + dLocalphiDot_dPhi + end subroutine source_damage_anisoDuctile_getRateAndItsTangent - module subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - integer, intent(in) :: & - phase, & !< phase ID of element - constituent !< position of element within its phase instance - real(pReal), intent(in) :: & - phi !< damage parameter - real(pReal), intent(out) :: & - localphiDot, & - dLocalphiDot_dPhi - end subroutine source_damage_isoBrittle_getRateAndItsTangent + module subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) + integer, intent(in) :: & + phase, & !< phase ID of element + constituent !< position of element within its phase instance + real(pReal), intent(in) :: & + phi !< damage parameter + real(pReal), intent(out) :: & + localphiDot, & + dLocalphiDot_dPhi + end subroutine source_damage_isoBrittle_getRateAndItsTangent - module subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - integer, intent(in) :: & - phase, & !< phase ID of element - constituent !< position of element within its phase instance - real(pReal), intent(in) :: & - phi !< damage parameter - real(pReal), intent(out) :: & - localphiDot, & - dLocalphiDot_dPhi - end subroutine source_damage_isoDuctile_getRateAndItsTangent + module subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) + integer, intent(in) :: & + phase, & !< phase ID of element + constituent !< position of element within its phase instance + real(pReal), intent(in) :: & + phi !< damage parameter + real(pReal), intent(out) :: & + localphiDot, & + dLocalphiDot_dPhi + end subroutine source_damage_isoDuctile_getRateAndItsTangent - module subroutine source_damage_anisoBrittle_results(phase,group) - integer, intent(in) :: phase - character(len=*), intent(in) :: group - end subroutine source_damage_anisoBrittle_results + module subroutine anisobrittle_results(phase,group) + integer, intent(in) :: phase + character(len=*), intent(in) :: group + end subroutine anisobrittle_results - module subroutine source_damage_anisoDuctile_results(phase,group) - integer, intent(in) :: phase - character(len=*), intent(in) :: group - end subroutine source_damage_anisoDuctile_results + module subroutine anisoductile_results(phase,group) + integer, intent(in) :: phase + character(len=*), intent(in) :: group + end subroutine anisoductile_results - module subroutine source_damage_isoBrittle_results(phase,group) - integer, intent(in) :: phase - character(len=*), intent(in) :: group - end subroutine source_damage_isoBrittle_results + module subroutine isobrittle_results(phase,group) + integer, intent(in) :: phase + character(len=*), intent(in) :: group + end subroutine isobrittle_results - module subroutine source_damage_isoDuctile_results(phase,group) - integer, intent(in) :: phase - character(len=*), intent(in) :: group - end subroutine source_damage_isoDuctile_results + module subroutine isoductile_results(phase,group) + integer, intent(in) :: phase + character(len=*), intent(in) :: group + end subroutine isoductile_results end interface @@ -164,19 +155,21 @@ module subroutine damage_init class(tNode), pointer :: & phases, & phase, & - sources, & - kinematics + sources + + + print'(/,a)', ' <<<+- phase:damage init -+>>>' phases => config_material%get('phase') allocate(current(phases%length)) allocate(damageState (phases%length)) - allocate(phase_Nsources(phases%length),source = 0) ! same for kinematics + allocate(phase_Nsources(phases%length),source = 0) do ph = 1,phases%length - Nconstituents = count(material_phaseAt == ph) * discretization_nIPs + Nconstituents = count(material_phaseAt2 == ph) allocate(current(ph)%phi(Nconstituents),source=1.0_pReal) allocate(current(ph)%d_phi_d_dot_phi(Nconstituents),source=0.0_pReal) @@ -191,26 +184,10 @@ module subroutine damage_init ! initialize source mechanisms if(maxval(phase_Nsources) /= 0) then - where(source_damage_isoBrittle_init (maxval(phase_Nsources))) phase_source = DAMAGE_ISOBRITTLE_ID - where(source_damage_isoDuctile_init (maxval(phase_Nsources))) phase_source = DAMAGE_ISODUCTILE_ID - where(source_damage_anisoBrittle_init (maxval(phase_Nsources))) phase_source = DAMAGE_ANISOBRITTLE_ID - where(source_damage_anisoDuctile_init (maxval(phase_Nsources))) phase_source = DAMAGE_ANISODUCTILE_ID - endif - -!-------------------------------------------------------------------------------------------------- -! initialize kinematic mechanisms - allocate(phase_Nkinematics(phases%length),source = 0) - do ph = 1,phases%length - phase => phases%get(ph) - kinematics => phase%get('kinematics',defaultVal=emptyList) - phase_Nkinematics(ph) = kinematics%length - enddo - - allocate(phase_kinematics(maxval(phase_Nkinematics),phases%length), source = KINEMATICS_undefined_ID) - - if(maxval(phase_Nkinematics) /= 0) then - where(kinematics_cleavage_opening_init(maxval(phase_Nkinematics))) phase_kinematics = KINEMATICS_cleavage_opening_ID - where(kinematics_slipplane_opening_init(maxval(phase_Nkinematics))) phase_kinematics = KINEMATICS_slipplane_opening_ID + where(isobrittle_init (maxval(phase_Nsources))) phase_source = DAMAGE_ISOBRITTLE_ID + where(isoductile_init (maxval(phase_Nsources))) phase_source = DAMAGE_ISODUCTILE_ID + where(anisobrittle_init (maxval(phase_Nsources))) phase_source = DAMAGE_ANISOBRITTLE_ID + where(anisoductile_init (maxval(phase_Nsources))) phase_source = DAMAGE_ANISODUCTILE_ID endif end subroutine damage_init @@ -234,30 +211,30 @@ module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi localphiDot, & dLocalphiDot_dPhi integer :: & - phase, & - grain, & - source, & - constituent + ph, & + co, & + so, & + me phiDot = 0.0_pReal dPhiDot_dPhi = 0.0_pReal - do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - phase = material_phaseAt(grain,el) - constituent = material_phasememberAt(grain,ip,el) - do source = 1, phase_Nsources(phase) - select case(phase_source(source,phase)) + do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) + ph = material_phaseAt(co,el) + me = material_phasememberAt(co,ip,el) + do so = 1, phase_Nsources(ph) + select case(phase_source(so,ph)) case (DAMAGE_ISOBRITTLE_ID) - call source_damage_isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) + call source_damage_isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, ph, me) case (DAMAGE_ISODUCTILE_ID) - call source_damage_isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) + call source_damage_isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, ph, me) case (DAMAGE_ANISOBRITTLE_ID) - call source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) + call source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) case (DAMAGE_ANISODUCTILE_ID) - call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) + call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) case default localphiDot = 0.0_pReal @@ -395,16 +372,16 @@ module subroutine damage_results(group,ph) sourceType: select case (phase_source(so,ph)) case (DAMAGE_ISOBRITTLE_ID) sourceType - call source_damage_isoBrittle_results(ph,group//'sources/') + call isobrittle_results(ph,group//'sources/') case (DAMAGE_ISODUCTILE_ID) sourceType - call source_damage_isoDuctile_results(ph,group//'sources/') + call isoductile_results(ph,group//'sources/') case (DAMAGE_ANISOBRITTLE_ID) sourceType - call source_damage_anisoBrittle_results(ph,group//'sources/') + call anisobrittle_results(ph,group//'sources/') case (DAMAGE_ANISODUCTILE_ID) sourceType - call source_damage_anisoDuctile_results(ph,group//'sources/') + call anisoductile_results(ph,group//'sources/') end select sourceType @@ -436,13 +413,13 @@ function constitutive_damage_collectDotState(co,ip,el,ph,me) result(broken) sourceType: select case (phase_source(so,ph)) case (DAMAGE_ISODUCTILE_ID) sourceType - call source_damage_isoDuctile_dotState(co, ip, el) + call isoductile_dotState(co, ip, el) case (DAMAGE_ANISODUCTILE_ID) sourceType - call source_damage_anisoDuctile_dotState(co, ip, el) + call anisoductile_dotState(co, ip, el) case (DAMAGE_ANISOBRITTLE_ID) sourceType - call source_damage_anisoBrittle_dotState(mech_S(ph,me),co, ip, el) ! correct stress? + call anisobrittle_dotState(mech_S(ph,me),co, ip, el) ! correct stress? end select sourceType @@ -542,7 +519,7 @@ end subroutine constitutive_damage_set_phi module function constitutive_damage_get_phi(co,ip,el) result(phi) - + integer, intent(in) :: co, ip, el real(pReal) :: phi @@ -551,4 +528,4 @@ module function constitutive_damage_get_phi(co,ip,el) result(phi) end function constitutive_damage_get_phi -end submodule constitutive_damage +end submodule damagee diff --git a/src/source_damage_anisoBrittle.f90 b/src/phase_damage_anisobrittle.f90 similarity index 94% rename from src/source_damage_anisoBrittle.f90 rename to src/phase_damage_anisobrittle.f90 index 8da37c5d5..2ae5ca951 100644 --- a/src/source_damage_anisoBrittle.f90 +++ b/src/phase_damage_anisobrittle.f90 @@ -4,7 +4,7 @@ !> @brief material subroutine incorporating anisotropic brittle damage source mechanism !> @details to be done !-------------------------------------------------------------------------------------------------- -submodule (constitutive:constitutive_damage) source_damage_anisoBrittle +submodule (phase:damagee) anisobrittle integer, dimension(:), allocatable :: & source_damage_anisoBrittle_offset, & !< which source is my current source mechanism? @@ -35,7 +35,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -module function source_damage_anisoBrittle_init(source_length) result(mySources) +module function anisobrittle_init(source_length) result(mySources) integer, intent(in) :: source_length logical, dimension(:,:), allocatable :: mySources @@ -49,7 +49,7 @@ module function source_damage_anisoBrittle_init(source_length) result(mySources) integer, dimension(:), allocatable :: N_cl character(len=pStringLen) :: extmsg = '' - print'(/,a)', ' <<<+- source_damage_anisoBrittle init -+>>>' + print'(/,a)', ' <<<+- phase:damage:anisobrittle init -+>>>' mySources = source_active('damage_anisoBrittle',source_length) Ninstances = count(mySources) @@ -114,13 +114,13 @@ module function source_damage_anisoBrittle_init(source_length) result(mySources) enddo enddo -end function source_damage_anisoBrittle_init +end function anisobrittle_init !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -module subroutine source_damage_anisoBrittle_dotState(S, co, ip, el) +module subroutine anisobrittle_dotState(S, co, ip, el) integer, intent(in) :: & co, & !< component-ID of integration point @@ -163,7 +163,7 @@ module subroutine source_damage_anisoBrittle_dotState(S, co, ip, el) enddo end associate -end subroutine source_damage_anisoBrittle_dotState +end subroutine anisobrittle_dotState !-------------------------------------------------------------------------------------------------- @@ -196,7 +196,7 @@ end subroutine source_damage_anisoBrittle_getRateAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- -module subroutine source_damage_anisoBrittle_results(phase,group) +module subroutine anisobrittle_results(phase,group) integer, intent(in) :: phase character(len=*), intent(in) :: group @@ -213,6 +213,6 @@ module subroutine source_damage_anisoBrittle_results(phase,group) enddo outputsLoop end associate -end subroutine source_damage_anisoBrittle_results +end subroutine anisobrittle_results -end submodule source_damage_anisoBrittle +end submodule anisobrittle diff --git a/src/source_damage_anisoDuctile.f90 b/src/phase_damage_anisoductile.f90 similarity index 93% rename from src/source_damage_anisoDuctile.f90 rename to src/phase_damage_anisoductile.f90 index b5e6fd17c..8d5904e6b 100644 --- a/src/source_damage_anisoDuctile.f90 +++ b/src/phase_damage_anisoductile.f90 @@ -4,7 +4,7 @@ !> @brief material subroutine incorporating anisotropic ductile damage source mechanism !> @details to be done !-------------------------------------------------------------------------------------------------- -submodule(constitutive:constitutive_damage) source_damage_anisoDuctile +submodule(phase:damagee) anisoductile integer, dimension(:), allocatable :: & source_damage_anisoDuctile_offset, & !< which source is my current damage mechanism? @@ -28,7 +28,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -module function source_damage_anisoDuctile_init(source_length) result(mySources) +module function anisoductile_init(source_length) result(mySources) integer, intent(in) :: source_length logical, dimension(:,:), allocatable :: mySources @@ -44,7 +44,7 @@ module function source_damage_anisoDuctile_init(source_length) result(mySources) integer, dimension(:), allocatable :: N_sl character(len=pStringLen) :: extmsg = '' - print'(/,a)', ' <<<+- source_damage_anisoDuctile init -+>>>' + print'(/,a)', ' <<<+- phase:damage:anisoductile init -+>>>' mySources = source_active('damage_anisoDuctile',source_length) Ninstances = count(mySources) @@ -101,13 +101,13 @@ module function source_damage_anisoDuctile_init(source_length) result(mySources) enddo -end function source_damage_anisoDuctile_init +end function anisoductile_init !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -module subroutine source_damage_anisoDuctile_dotState(co, ip, el) +module subroutine anisoductile_dotState(co, ip, el) integer, intent(in) :: & co, & !< component-ID of integration point @@ -132,7 +132,7 @@ module subroutine source_damage_anisoDuctile_dotState(co, ip, el) = sum(plasticState(ph)%slipRate(:,me)/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit) end associate -end subroutine source_damage_anisoDuctile_dotState +end subroutine anisoductile_dotState !-------------------------------------------------------------------------------------------------- @@ -165,7 +165,7 @@ end subroutine source_damage_anisoDuctile_getRateAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- -module subroutine source_damage_anisoDuctile_results(phase,group) +module subroutine anisoductile_results(phase,group) integer, intent(in) :: phase character(len=*), intent(in) :: group @@ -182,6 +182,6 @@ module subroutine source_damage_anisoDuctile_results(phase,group) enddo outputsLoop end associate -end subroutine source_damage_anisoDuctile_results +end subroutine anisoductile_results -end submodule source_damage_anisoDuctile +end submodule anisoductile diff --git a/src/source_damage_isoBrittle.f90 b/src/phase_damage_isobrittle.f90 similarity index 94% rename from src/source_damage_isoBrittle.f90 rename to src/phase_damage_isobrittle.f90 index cff993e7f..091377171 100644 --- a/src/source_damage_isoBrittle.f90 +++ b/src/phase_damage_isobrittle.f90 @@ -4,7 +4,7 @@ !> @brief material subroutine incoprorating isotropic brittle damage source mechanism !> @details to be done !-------------------------------------------------------------------------------------------------- -submodule(constitutive:constitutive_damage) source_damage_isoBrittle +submodule(phase:damagee) isobrittle integer, dimension(:), allocatable :: & source_damage_isoBrittle_offset, & @@ -26,7 +26,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -module function source_damage_isoBrittle_init(source_length) result(mySources) +module function isobrittle_init(source_length) result(mySources) integer, intent(in) :: source_length logical, dimension(:,:), allocatable :: mySources @@ -39,7 +39,7 @@ module function source_damage_isoBrittle_init(source_length) result(mySources) integer :: Ninstances,sourceOffset,Nconstituents,p character(len=pStringLen) :: extmsg = '' - print'(/,a)', ' <<<+- source_damage_isoBrittle init -+>>>' + print'(/,a)', ' <<<+- phase:damage:isobrittle init -+>>>' mySources = source_active('damage_isoBrittle',source_length) Ninstances = count(mySources) @@ -88,7 +88,7 @@ module function source_damage_isoBrittle_init(source_length) result(mySources) enddo -end function source_damage_isoBrittle_init +end function isobrittle_init !-------------------------------------------------------------------------------------------------- @@ -161,7 +161,7 @@ end subroutine source_damage_isoBrittle_getRateAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- -module subroutine source_damage_isoBrittle_results(phase,group) +module subroutine isobrittle_results(phase,group) integer, intent(in) :: phase character(len=*), intent(in) :: group @@ -178,6 +178,6 @@ module subroutine source_damage_isoBrittle_results(phase,group) enddo outputsLoop end associate -end subroutine source_damage_isoBrittle_results +end subroutine isobrittle_results -end submodule source_damage_isoBrittle +end submodule isobrittle diff --git a/src/source_damage_isoDuctile.f90 b/src/phase_damage_isoductile.f90 similarity index 92% rename from src/source_damage_isoDuctile.f90 rename to src/phase_damage_isoductile.f90 index 3cebc474a..6f0a2d0fb 100644 --- a/src/source_damage_isoDuctile.f90 +++ b/src/phase_damage_isoductile.f90 @@ -4,7 +4,7 @@ !> @brief material subroutine incorporating isotropic ductile damage source mechanism !> @details to be done !-------------------------------------------------------------------------------------------------- -submodule (constitutive:constitutive_damage) source_damage_isoDuctile +submodule(phase:damagee) isoductile integer, dimension(:), allocatable :: & source_damage_isoDuctile_offset, & !< which source is my current damage mechanism? @@ -28,7 +28,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -module function source_damage_isoDuctile_init(source_length) result(mySources) +module function isoductile_init(source_length) result(mySources) integer, intent(in) :: source_length logical, dimension(:,:), allocatable :: mySources @@ -41,7 +41,7 @@ module function source_damage_isoDuctile_init(source_length) result(mySources) integer :: Ninstances,sourceOffset,Nconstituents,p character(len=pStringLen) :: extmsg = '' - print'(/,a)', ' <<<+- source_damage_isoDuctile init -+>>>' + print'(/,a)', ' <<<+- phase:damage:isoductile init -+>>>' mySources = source_active('damage_isoDuctile',source_length) Ninstances = count(mySources) @@ -92,13 +92,13 @@ module function source_damage_isoDuctile_init(source_length) result(mySources) enddo -end function source_damage_isoDuctile_init +end function isoductile_init !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -module subroutine source_damage_isoDuctile_dotState(co, ip, el) +module subroutine isoductile_dotState(co, ip, el) integer, intent(in) :: & co, & !< component-ID of integration point @@ -123,7 +123,7 @@ module subroutine source_damage_isoDuctile_dotState(co, ip, el) sum(plasticState(ph)%slipRate(:,me))/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit end associate -end subroutine source_damage_isoDuctile_dotState +end subroutine isoductile_dotState !-------------------------------------------------------------------------------------------------- @@ -156,7 +156,7 @@ end subroutine source_damage_isoDuctile_getRateAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- -module subroutine source_damage_isoDuctile_results(phase,group) +module subroutine isoductile_results(phase,group) integer, intent(in) :: phase character(len=*), intent(in) :: group @@ -173,6 +173,6 @@ module subroutine source_damage_isoDuctile_results(phase,group) enddo outputsLoop end associate -end subroutine source_damage_isoDuctile_results +end subroutine isoductile_results -end submodule source_damage_isoDuctile +end submodule isoductile diff --git a/src/constitutive_mech.f90 b/src/phase_mechanics.f90 similarity index 71% rename from src/constitutive_mech.f90 rename to src/phase_mechanics.f90 index be933ef4f..a84c1a385 100644 --- a/src/constitutive_mech.f90 +++ b/src/phase_mechanics.f90 @@ -1,7 +1,7 @@ !---------------------------------------------------------------------------------------------------- !> @brief internal microstructure state for all plasticity constitutive models !---------------------------------------------------------------------------------------------------- -submodule(constitutive) constitutive_mech +submodule(phase) mechanics enum, bind(c); enumerator :: & ELASTICITY_UNDEFINED_ID, & @@ -15,9 +15,15 @@ submodule(constitutive) constitutive_mech PLASTICITY_KINEHARDENING_ID, & PLASTICITY_DISLOTWIN_ID, & PLASTICITY_DISLOTUNGSTEN_ID, & - PLASTICITY_NONLOCAL_ID + PLASTICITY_NONLOCAL_ID, & + KINEMATICS_UNDEFINED_ID, & + KINEMATICS_CLEAVAGE_OPENING_ID, & + KINEMATICS_SLIPPLANE_OPENING_ID, & + KINEMATICS_THERMAL_EXPANSION_ID end enum + integer(kind(KINEMATICS_UNDEFINED_ID)), dimension(:,:), allocatable :: & + phase_kinematics integer(kind(ELASTICITY_UNDEFINED_ID)), dimension(:), allocatable :: & phase_elasticity !< elasticity of each phase integer(kind(STIFFNESS_DEGRADATION_UNDEFINED_ID)), dimension(:,:), allocatable :: & @@ -48,225 +54,85 @@ submodule(constitutive) constitutive_mech interface - module function plastic_none_init() result(myPlasticity) - logical, dimension(:), allocatable :: & - myPlasticity - end function plastic_none_init + module subroutine eigendeformation_init(phases) + class(tNode), pointer :: phases + end subroutine eigendeformation_init - module function plastic_isotropic_init() result(myPlasticity) - logical, dimension(:), allocatable :: & - myPlasticity - end function plastic_isotropic_init + module subroutine plastic_init + end subroutine plastic_init - module function plastic_phenopowerlaw_init() result(myPlasticity) - logical, dimension(:), allocatable :: & - myPlasticity - end function plastic_phenopowerlaw_init - - module function plastic_kinehardening_init() result(myPlasticity) - logical, dimension(:), allocatable :: & - myPlasticity - end function plastic_kinehardening_init - - module function plastic_dislotwin_init() result(myPlasticity) - logical, dimension(:), allocatable :: & - myPlasticity - end function plastic_dislotwin_init - - module function plastic_dislotungsten_init() result(myPlasticity) - logical, dimension(:), allocatable :: & - myPlasticity - end function plastic_dislotungsten_init - - module function plastic_nonlocal_init() result(myPlasticity) - logical, dimension(:), allocatable :: & - myPlasticity - end function plastic_nonlocal_init - - - module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) + module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,me) real(pReal), dimension(3,3), intent(out) :: & - Lp !< plastic velocity gradient - real(pReal), dimension(3,3,3,3), intent(out) :: & - dLp_dMp !< derivative of Lp with respect to the Mandel stress - + Li !< inleastic velocity gradient + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLi_dMi !< derivative of Li with respect to Mandel stress real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress + Mi !< Mandel stress integer, intent(in) :: & instance, & - of - end subroutine plastic_isotropic_LpAndItsTangent + me + end subroutine plastic_isotropic_LiAndItsTangent - pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) - real(pReal), dimension(3,3), intent(out) :: & - Lp !< plastic velocity gradient - real(pReal), dimension(3,3,3,3), intent(out) :: & - dLp_dMp !< derivative of Lp with respect to the Mandel stress - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - end subroutine plastic_phenopowerlaw_LpAndItsTangent + module function plastic_dotState(subdt,co,ip,el,ph,me) result(broken) - pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) - real(pReal), dimension(3,3), intent(out) :: & - Lp !< plastic velocity gradient - real(pReal), dimension(3,3,3,3), intent(out) :: & - dLp_dMp !< derivative of Lp with respect to the Mandel stress - - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - end subroutine plastic_kinehardening_LpAndItsTangent - - module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) - real(pReal), dimension(3,3), intent(out) :: & - Lp !< plastic velocity gradient - real(pReal), dimension(3,3,3,3), intent(out) :: & - dLp_dMp !< derivative of Lp with respect to the Mandel stress - - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - real(pReal), intent(in) :: & - T - integer, intent(in) :: & - instance, & - of - end subroutine plastic_dislotwin_LpAndItsTangent - - pure module subroutine plastic_dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) - real(pReal), dimension(3,3), intent(out) :: & - Lp !< plastic velocity gradient - real(pReal), dimension(3,3,3,3), intent(out) :: & - dLp_dMp !< derivative of Lp with respect to the Mandel stress - - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - real(pReal), intent(in) :: & - T - integer, intent(in) :: & - instance, & - of - end subroutine plastic_dislotungsten_LpAndItsTangent - - module subroutine plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp, & - Mp,Temperature,instance,of,ip,el) - real(pReal), dimension(3,3), intent(out) :: & - Lp !< plastic velocity gradient - real(pReal), dimension(3,3,3,3), intent(out) :: & - dLp_dMp !< derivative of Lp with respect to the Mandel stress - - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - real(pReal), intent(in) :: & - Temperature - integer, intent(in) :: & - instance, & - of, & - ip, & !< current integration point - el !< current element number - end subroutine plastic_nonlocal_LpAndItsTangent - - module subroutine plastic_isotropic_dotState(Mp,instance,of) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - end subroutine plastic_isotropic_dotState - - module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - end subroutine plastic_phenopowerlaw_dotState - - module subroutine plastic_kinehardening_dotState(Mp,instance,of) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - end subroutine plastic_kinehardening_dotState - - module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - real(pReal), intent(in) :: & - T - integer, intent(in) :: & - instance, & - of - end subroutine plastic_dislotwin_dotState - - module subroutine plastic_disloTungsten_dotState(Mp,T,instance,of) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - real(pReal), intent(in) :: & - T - integer, intent(in) :: & - instance, & - of - end subroutine plastic_disloTungsten_dotState - - module subroutine plastic_nonlocal_dotState(Mp,Temperature,timestep,instance,of,ip,el) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< MandelStress - real(pReal), intent(in) :: & - Temperature, & !< temperature - timestep !< substepped crystallite time increment integer, intent(in) :: & - instance, & - of, & - ip, & !< current integration point - el !< current element number - end subroutine plastic_nonlocal_dotState + co, & !< component-ID of integration point + ip, & !< integration point + el, & !< element + ph, & + me + real(pReal), intent(in) :: & + subdt !< timestep + logical :: broken + end function plastic_dotState - - module subroutine plastic_dislotwin_dependentState(T,instance,of) - integer, intent(in) :: & - instance, & - of - real(pReal), intent(in) :: & - T - end subroutine plastic_dislotwin_dependentState - - module subroutine plastic_dislotungsten_dependentState(instance,of) - integer, intent(in) :: & - instance, & - of - end subroutine plastic_dislotungsten_dependentState - - module subroutine plastic_nonlocal_dependentState(instance, of, ip, el) + module function plastic_deltaState(co, ip, el, ph, me) result(broken) integer, intent(in) :: & - instance, & - of, & - ip, & !< current integration point - el !< current element number - end subroutine plastic_nonlocal_dependentState + co, & !< component-ID of integration point + ip, & !< integration point + el, & !< element + ph, & + me + logical :: & + broken + end function plastic_deltaState - module subroutine plastic_kinehardening_deltaState(Mp,instance,of) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - end subroutine plastic_kinehardening_deltaState + module subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & + S, Fi, co, ip, el) - module subroutine plastic_nonlocal_deltaState(Mp,instance,of,ip,el) - real(pReal), dimension(3,3), intent(in) :: & - Mp integer, intent(in) :: & - instance, & - of, & - ip, & - el - end subroutine plastic_nonlocal_deltaState + co, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + S !< 2nd Piola-Kirchhoff stress + real(pReal), intent(in), dimension(3,3) :: & + Fi !< intermediate deformation gradient + real(pReal), intent(out), dimension(3,3) :: & + Li !< intermediate velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLi_dS, & !< derivative of Li with respect to S + dLi_dFi + + end subroutine constitutive_LiAndItsTangents + + + module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & + S, Fi, co, ip, el) + integer, intent(in) :: & + co, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + S, & !< 2nd Piola-Kirchhoff stress + Fi !< intermediate deformation gradient + real(pReal), intent(out), dimension(3,3) :: & + Lp !< plastic velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLp_dS, & + dLp_dFi !< derivative of Lp with respect to Fi + end subroutine plastic_LpAndItsTangents + module subroutine plastic_isotropic_results(instance,group) integer, intent(in) :: instance @@ -340,7 +206,7 @@ module subroutine mech_init(phases) elastic, & stiffDegradation - print'(/,a)', ' <<<+- constitutive_mech init -+>>>' + print'(/,a)', ' <<<+- phase:mechanics init -+>>>' !------------------------------------------------------------------------------------------------- ! initialize elasticity (hooke) !ToDO: Maybe move to elastic submodule along with function homogenizedC? @@ -444,13 +310,7 @@ module subroutine mech_init(phases) allocate(phase_plasticityInstance(phases%length),source = 0) allocate(phase_localPlasticity(phases%length), source=.true.) - where(plastic_none_init()) phase_plasticity = PLASTICITY_NONE_ID - where(plastic_isotropic_init()) phase_plasticity = PLASTICITY_ISOTROPIC_ID - where(plastic_phenopowerlaw_init()) phase_plasticity = PLASTICITY_PHENOPOWERLAW_ID - where(plastic_kinehardening_init()) phase_plasticity = PLASTICITY_KINEHARDENING_ID - where(plastic_dislotwin_init()) phase_plasticity = PLASTICITY_DISLOTWIN_ID - where(plastic_dislotungsten_init()) phase_plasticity = PLASTICITY_DISLOTUNGSTEN_ID - where(plastic_nonlocal_init()) phase_plasticity = PLASTICITY_NONLOCAL_ID + call plastic_init() do ph = 1, phases%length phase_elasticityInstance(ph) = count(phase_elasticity(1:ph) == phase_elasticity(ph)) @@ -481,36 +341,13 @@ module subroutine mech_init(phases) end select + + call eigendeformation_init(phases) + + end subroutine mech_init -!-------------------------------------------------------------------------------------------------- -!> @brief checks if a plastic module is active or not -!-------------------------------------------------------------------------------------------------- -function plastic_active(plastic_label) result(active_plastic) - - character(len=*), intent(in) :: plastic_label !< type of plasticity model - logical, dimension(:), allocatable :: active_plastic - - class(tNode), pointer :: & - phases, & - phase, & - mech, & - pl - integer :: ph - - phases => config_material%get('phase') - allocate(active_plastic(phases%length), source = .false. ) - do ph = 1, phases%length - phase => phases%get(ph) - mech => phase%get('mechanics') - pl => mech%get('plasticity') - if(pl%get_asString('type') == plastic_label) active_plastic(ph) = .true. - enddo - -end function plastic_active - - !-------------------------------------------------------------------------------------------------- !> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to !> the elastic and intermediate deformation gradients using Hooke's law @@ -559,216 +396,6 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & end subroutine constitutive_hooke_SandItsTangents -!-------------------------------------------------------------------------------------------------- -!> @brief calls microstructure function of the different plasticity constitutive models -!-------------------------------------------------------------------------------------------------- -module subroutine constitutive_plastic_dependentState(co, ip, el) - - integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element - - integer :: & - ph, & - instance, me - - - ph = material_phaseAt(co,el) - me = material_phasememberAt(co,ip,el) - instance = phase_plasticityInstance(ph) - - plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) - - case (PLASTICITY_DISLOTWIN_ID) plasticityType - call plastic_dislotwin_dependentState(thermal_T(ph,me),instance,me) - - case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType - call plastic_dislotungsten_dependentState(instance,me) - - case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_dependentState(instance,me,ip,el) - - end select plasticityType - -end subroutine constitutive_plastic_dependentState - - -!-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the velocity gradient -! ToDo: Discuss whether it makes sense if crystallite handles the configuration conversion, i.e. -! Mp in, dLp_dMp out -!-------------------------------------------------------------------------------------------------- -subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & - S, Fi, co, ip, el) - integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element - real(pReal), intent(in), dimension(3,3) :: & - S, & !< 2nd Piola-Kirchhoff stress - Fi !< intermediate deformation gradient - real(pReal), intent(out), dimension(3,3) :: & - Lp !< plastic velocity gradient - real(pReal), intent(out), dimension(3,3,3,3) :: & - dLp_dS, & - dLp_dFi !< derivative of Lp with respect to Fi - - real(pReal), dimension(3,3,3,3) :: & - dLp_dMp !< derivative of Lp with respect to Mandel stress - real(pReal), dimension(3,3) :: & - Mp !< Mandel stress work conjugate with Lp - integer :: & - i, j, instance, me, ph - - - Mp = matmul(matmul(transpose(Fi),Fi),S) - me = material_phasememberAt(co,ip,el) - ph = material_phaseAt(co,el) - instance = phase_plasticityInstance(ph) - - plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) - - case (PLASTICITY_NONE_ID) plasticityType - Lp = 0.0_pReal - dLp_dMp = 0.0_pReal - - case (PLASTICITY_ISOTROPIC_ID) plasticityType - call plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me) - - case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType - call plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me) - - case (PLASTICITY_KINEHARDENING_ID) plasticityType - call plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me) - - case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),instance,me,ip,el) - - case (PLASTICITY_DISLOTWIN_ID) plasticityType - call plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),instance,me) - - case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType - call plastic_dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),instance,me) - - end select plasticityType - - do i=1,3; do j=1,3 - dLp_dFi(i,j,1:3,1:3) = matmul(matmul(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + & - matmul(matmul(Fi,dLp_dMp(i,j,1:3,1:3)),S) - dLp_dS(i,j,1:3,1:3) = matmul(matmul(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi) - enddo; enddo - -end subroutine constitutive_plastic_LpAndItsTangents - - -!-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the rate of change of microstructure -!-------------------------------------------------------------------------------------------------- -function mech_collectDotState(subdt,co,ip,el,ph,me) result(broken) - - integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el, & !< element - ph, & - me - real(pReal), intent(in) :: & - subdt !< timestep - real(pReal), dimension(3,3) :: & - Mp - integer :: & - instance - logical :: broken - - - instance = phase_plasticityInstance(ph) - - Mp = matmul(matmul(transpose(constitutive_mech_Fi(ph)%data(1:3,1:3,me)),& - constitutive_mech_Fi(ph)%data(1:3,1:3,me)),constitutive_mech_S(ph)%data(1:3,1:3,me)) - - plasticityType: select case (phase_plasticity(ph)) - - case (PLASTICITY_ISOTROPIC_ID) plasticityType - call plastic_isotropic_dotState(Mp,instance,me) - - case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType - call plastic_phenopowerlaw_dotState(Mp,instance,me) - - case (PLASTICITY_KINEHARDENING_ID) plasticityType - call plastic_kinehardening_dotState(Mp,instance,me) - - case (PLASTICITY_DISLOTWIN_ID) plasticityType - call plastic_dislotwin_dotState(Mp,thermal_T(ph,me),instance,me) - - case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType - call plastic_disloTungsten_dotState(Mp,thermal_T(ph,me),instance,me) - - case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_dotState(Mp,thermal_T(ph,me),subdt,instance,me,ip,el) - end select plasticityType - broken = any(IEEE_is_NaN(plasticState(ph)%dotState(:,me))) - - -end function mech_collectDotState - - -!-------------------------------------------------------------------------------------------------- -!> @brief for constitutive models having an instantaneous change of state -!> will return false if delta state is not needed/supported by the constitutive model -!-------------------------------------------------------------------------------------------------- -function constitutive_deltaState(co, ip, el, ph, of) result(broken) - - integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el, & !< element - ph, & - of - logical :: & - broken - - real(pReal), dimension(3,3) :: & - Mp - integer :: & - instance, & - myOffset, & - mySize - - - Mp = matmul(matmul(transpose(constitutive_mech_Fi(ph)%data(1:3,1:3,of)),& - constitutive_mech_Fi(ph)%data(1:3,1:3,of)),constitutive_mech_S(ph)%data(1:3,1:3,of)) - instance = phase_plasticityInstance(ph) - - plasticityType: select case (phase_plasticity(ph)) - - case (PLASTICITY_KINEHARDENING_ID) plasticityType - call plastic_kinehardening_deltaState(Mp,instance,of) - broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,of))) - - case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_deltaState(Mp,instance,of,ip,el) - broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,of))) - - case default - broken = .false. - - end select plasticityType - - if(.not. broken) then - select case(phase_plasticity(ph)) - case (PLASTICITY_NONLOCAL_ID,PLASTICITY_KINEHARDENING_ID) - - myOffset = plasticState(ph)%offsetDeltaState - mySize = plasticState(ph)%sizeDeltaState - plasticState(ph)%state(myOffset + 1:myOffset + mySize,of) = & - plasticState(ph)%state(myOffset + 1:myOffset + mySize,of) + plasticState(ph)%deltaState(1:mySize,of) - end select - endif - -end function constitutive_deltaState - - module subroutine mech_results(group,ph) character(len=*), intent(in) :: group @@ -862,7 +489,6 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken) ierr, & ! error indicator for LAPACK o, & p, & - m, & ph, & me, & jacoCounterLp, & @@ -875,7 +501,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken) ph = material_phaseAt(co,el) me = material_phaseMemberAt(co,ip,el) - call constitutive_plastic_dependentState(co,ip,el) + call plastic_dependentState(co,ip,el) Lpguess = constitutive_mech_Lp(ph)%data(1:3,1:3,me) ! take as first guess Liguess = constitutive_mech_Li(ph)%data(1:3,1:3,me) ! take as first guess @@ -915,7 +541,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken) call constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & Fe, Fi_new, co, ip, el) - call constitutive_plastic_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, & + call plastic_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, & S, Fi_new, co, ip, el) !* update current residuum and check for convergence of loop @@ -1051,7 +677,7 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul ph = material_phaseAt(co,el) me = material_phaseMemberAt(co,ip,el) - broken = mech_collectDotState(Delta_t, co,ip,el,ph,me) + broken = plastic_dotState(Delta_t, co,ip,el,ph,me) if(broken) return sizeDotState = plasticState(ph)%sizeDotState @@ -1067,7 +693,7 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) if(broken) exit iteration - broken = mech_collectDotState(Delta_t, co,ip,el,ph,me) + broken = plastic_dotState(Delta_t, co,ip,el,ph,me) if(broken) exit iteration zeta = damper(plasticState(ph)%dotState(:,me),dotState(1:sizeDotState,1),& @@ -1080,7 +706,7 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul plasticState(ph)%state(1:sizeDotState,me) = plasticState(ph)%state(1:sizeDotState,me) & - r(1:sizeDotState) if (converged(r(1:sizeDotState),plasticState(ph)%state(1:sizeDotState,me),plasticState(ph)%atol(1:sizeDotState))) then - broken = constitutive_deltaState(co,ip,el,ph,me) + broken = plastic_deltaState(co,ip,el,ph,me) exit iteration endif @@ -1136,14 +762,14 @@ function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) res ph = material_phaseAt(co,el) me = material_phaseMemberAt(co,ip,el) - broken = mech_collectDotState(Delta_t, co,ip,el,ph,me) + broken = plastic_dotState(Delta_t, co,ip,el,ph,me) if(broken) return sizeDotState = plasticState(ph)%sizeDotState plasticState(ph)%state(1:sizeDotState,me) = subState0 & + plasticState(ph)%dotState(1:sizeDotState,me) * Delta_t - broken = constitutive_deltaState(co,ip,el,ph,me) + broken = plastic_deltaState(co,ip,el,ph,me) if(broken) return broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) @@ -1176,7 +802,7 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip ph = material_phaseAt(co,el) me = material_phaseMemberAt(co,ip,el) - broken = mech_collectDotState(Delta_t, co,ip,el,ph,me) + broken = plastic_dotState(Delta_t, co,ip,el,ph,me) if(broken) return sizeDotState = plasticState(ph)%sizeDotState @@ -1185,13 +811,13 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip plasticState(ph)%state(1:sizeDotState,me) = subState0 & + plasticState(ph)%dotstate(1:sizeDotState,me) * Delta_t - broken = constitutive_deltaState(co,ip,el,ph,me) + broken = plastic_deltaState(co,ip,el,ph,me) if(broken) return broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) if(broken) return - broken = mech_collectDotState(Delta_t, co,ip,el,ph,me) + broken = plastic_dotState(Delta_t, co,ip,el,ph,me) if(broken) return broken = .not. converged(residuum_plastic(1:sizeDotState) + 0.5_pReal * plasticState(ph)%dotState(:,me) * Delta_t, & @@ -1294,7 +920,7 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D ph = material_phaseAt(co,el) me = material_phaseMemberAt(co,ip,el) - broken = mech_collectDotState(Delta_t,co,ip,el,ph,me) + broken = plastic_dotState(Delta_t,co,ip,el,ph,me) if(broken) return sizeDotState = plasticState(ph)%sizeDotState @@ -1315,7 +941,7 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D broken = integrateStress(F_0 + (F - F_0) * Delta_t * C(stage),subFp0,subFi0,Delta_t * C(stage),co,ip,el) if(broken) exit - broken = mech_collectDotState(Delta_t*C(stage),co,ip,el,ph,me) + broken = plastic_dotState(Delta_t*C(stage),co,ip,el,ph,me) if(broken) exit enddo @@ -1334,7 +960,7 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D if(broken) return - broken = constitutive_deltaState(co,ip,el,ph,me) + broken = plastic_deltaState(co,ip,el,ph,me) if(broken) return broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) @@ -1351,7 +977,6 @@ subroutine crystallite_results(group,ph) integer, intent(in) :: ph integer :: ou - real(pReal), allocatable, dimension(:,:,:) :: selected_tensors real(pReal), allocatable, dimension(:,:) :: selected_rotations character(len=:), allocatable :: structureLabel @@ -1708,7 +1333,7 @@ module function constitutive_mech_dPdF(dt,co,ip,el) result(dPdF) dLidS = math_mul3333xx3333(dLidFi,dFidS) + dLidS endif - call constitutive_plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, & + call plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, & constitutive_mech_S(ph)%data(1:3,1:3,me), & constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el) dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS @@ -1880,86 +1505,4 @@ module subroutine constitutive_mech_setF(F,co,ip,el) end subroutine constitutive_mech_setF -!-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the velocity gradient -! ToDo: MD: S is Mi? -!-------------------------------------------------------------------------------------------------- -subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & - S, Fi, co, ip, el) - - integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element - real(pReal), intent(in), dimension(3,3) :: & - S !< 2nd Piola-Kirchhoff stress - real(pReal), intent(in), dimension(3,3) :: & - Fi !< intermediate deformation gradient - real(pReal), intent(out), dimension(3,3) :: & - Li !< intermediate velocity gradient - real(pReal), intent(out), dimension(3,3,3,3) :: & - dLi_dS, & !< derivative of Li with respect to S - dLi_dFi - - real(pReal), dimension(3,3) :: & - my_Li, & !< intermediate velocity gradient - FiInv, & - temp_33 - real(pReal), dimension(3,3,3,3) :: & - my_dLi_dS - real(pReal) :: & - detFi - integer :: & - k, i, j, & - instance, of, me, ph - - Li = 0.0_pReal - dLi_dS = 0.0_pReal - dLi_dFi = 0.0_pReal - - plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) - case (PLASTICITY_isotropic_ID) plasticityType - of = material_phasememberAt(co,ip,el) - instance = phase_plasticityInstance(material_phaseAt(co,el)) - call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,instance,of) - case default plasticityType - my_Li = 0.0_pReal - my_dLi_dS = 0.0_pReal - end select plasticityType - - Li = Li + my_Li - dLi_dS = dLi_dS + my_dLi_dS - - KinematicsLoop: do k = 1, phase_Nkinematics(material_phaseAt(co,el)) - kinematicsType: select case (phase_kinematics(k,material_phaseAt(co,el))) - case (KINEMATICS_cleavage_opening_ID) kinematicsType - call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el) - case (KINEMATICS_slipplane_opening_ID) kinematicsType - call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el) - case (KINEMATICS_thermal_expansion_ID) kinematicsType - me = material_phaseMemberAt(co,ip,el) - ph = material_phaseAt(co,el) - call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dS, ph,me) - case default kinematicsType - my_Li = 0.0_pReal - my_dLi_dS = 0.0_pReal - end select kinematicsType - Li = Li + my_Li - dLi_dS = dLi_dS + my_dLi_dS - enddo KinematicsLoop - - FiInv = math_inv33(Fi) - detFi = math_det33(Fi) - Li = matmul(matmul(Fi,Li),FiInv)*detFi !< push forward to intermediate configuration - temp_33 = matmul(FiInv,Li) - - do i = 1,3; do j = 1,3 - dLi_dS(1:3,1:3,i,j) = matmul(matmul(Fi,dLi_dS(1:3,1:3,i,j)),FiInv)*detFi - dLi_dFi(1:3,1:3,i,j) = dLi_dFi(1:3,1:3,i,j) + Li*FiInv(j,i) - dLi_dFi(1:3,i,1:3,j) = dLi_dFi(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i) - enddo; enddo - -end subroutine constitutive_LiAndItsTangents - -end submodule constitutive_mech - +end submodule mechanics diff --git a/src/phase_mechanics_eigendeformation.f90 b/src/phase_mechanics_eigendeformation.f90 new file mode 100644 index 000000000..45cfd82d3 --- /dev/null +++ b/src/phase_mechanics_eigendeformation.f90 @@ -0,0 +1,207 @@ +submodule(phase:mechanics) eigendeformation + + interface + module function kinematics_cleavage_opening_init(kinematics_length) result(myKinematics) + integer, intent(in) :: kinematics_length + logical, dimension(:,:), allocatable :: myKinematics + end function kinematics_cleavage_opening_init + + module function kinematics_slipplane_opening_init(kinematics_length) result(myKinematics) + integer, intent(in) :: kinematics_length + logical, dimension(:,:), allocatable :: myKinematics + end function kinematics_slipplane_opening_init + + module function kinematics_thermal_expansion_init(kinematics_length) result(myKinematics) + integer, intent(in) :: kinematics_length + logical, dimension(:,:), allocatable :: myKinematics + end function kinematics_thermal_expansion_init + + module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el) + integer, intent(in) :: & + co, & !< grain number + ip, & !< integration point number + el !< element number + real(pReal), intent(in), dimension(3,3) :: & + S + real(pReal), intent(out), dimension(3,3) :: & + Ld !< damage velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) + end subroutine kinematics_cleavage_opening_LiAndItsTangent + + module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el) + integer, intent(in) :: & + co, & !< grain number + ip, & !< integration point number + el !< element number + real(pReal), intent(in), dimension(3,3) :: & + S + real(pReal), intent(out), dimension(3,3) :: & + Ld !< damage velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) + end subroutine kinematics_slipplane_opening_LiAndItsTangent + + module subroutine thermalexpansion_LiAndItsTangent(Li, dLi_dTstar, ph,me) + integer, intent(in) :: ph, me + !< element number + real(pReal), intent(out), dimension(3,3) :: & + Li !< thermal velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero) + end subroutine thermalexpansion_LiAndItsTangent + + end interface + + +contains + + +module subroutine eigendeformation_init(phases) + + class(tNode), pointer :: & + phases + + integer :: & + ph + class(tNode), pointer :: & + phase, & + kinematics + + print'(/,a)', ' <<<+- phase:mechanics:eigendeformation init -+>>>' + +!-------------------------------------------------------------------------------------------------- +! initialize kinematic mechanisms + allocate(phase_Nkinematics(phases%length),source = 0) + do ph = 1,phases%length + phase => phases%get(ph) + kinematics => phase%get('kinematics',defaultVal=emptyList) + phase_Nkinematics(ph) = kinematics%length + enddo + + allocate(phase_kinematics(maxval(phase_Nkinematics),phases%length), source = KINEMATICS_undefined_ID) + + if(maxval(phase_Nkinematics) /= 0) then + where(kinematics_cleavage_opening_init(maxval(phase_Nkinematics))) phase_kinematics = KINEMATICS_cleavage_opening_ID + where(kinematics_slipplane_opening_init(maxval(phase_Nkinematics))) phase_kinematics = KINEMATICS_slipplane_opening_ID + where(kinematics_thermal_expansion_init(maxval(phase_Nkinematics))) phase_kinematics = KINEMATICS_thermal_expansion_ID + endif + +end subroutine eigendeformation_init + + + +!-------------------------------------------------------------------------------------------------- +!> @brief checks if a kinematic mechanism is active or not +!-------------------------------------------------------------------------------------------------- +function kinematics_active(kinematics_label,kinematics_length) result(active_kinematics) + + character(len=*), intent(in) :: kinematics_label !< name of kinematic mechanism + integer, intent(in) :: kinematics_length !< max. number of kinematics in system + logical, dimension(:,:), allocatable :: active_kinematics + + class(tNode), pointer :: & + phases, & + phase, & + kinematics, & + kinematics_type + integer :: p,k + + phases => config_material%get('phase') + allocate(active_kinematics(kinematics_length,phases%length), source = .false. ) + do p = 1, phases%length + phase => phases%get(p) + kinematics => phase%get('kinematics',defaultVal=emptyList) + do k = 1, kinematics%length + kinematics_type => kinematics%get(k) + if(kinematics_type%get_asString('type') == kinematics_label) active_kinematics(k,p) = .true. + enddo + enddo + + +end function kinematics_active + + +!-------------------------------------------------------------------------------------------------- +!> @brief contains the constitutive equation for calculating the velocity gradient +! ToDo: MD: S is Mi? +!-------------------------------------------------------------------------------------------------- +module subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & + S, Fi, co, ip, el) + + integer, intent(in) :: & + co, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + S !< 2nd Piola-Kirchhoff stress + real(pReal), intent(in), dimension(3,3) :: & + Fi !< intermediate deformation gradient + real(pReal), intent(out), dimension(3,3) :: & + Li !< intermediate velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLi_dS, & !< derivative of Li with respect to S + dLi_dFi + + real(pReal), dimension(3,3) :: & + my_Li, & !< intermediate velocity gradient + FiInv, & + temp_33 + real(pReal), dimension(3,3,3,3) :: & + my_dLi_dS + real(pReal) :: & + detFi + integer :: & + k, i, j, & + instance, of, me, ph + + Li = 0.0_pReal + dLi_dS = 0.0_pReal + dLi_dFi = 0.0_pReal + + plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) + case (PLASTICITY_isotropic_ID) plasticityType + of = material_phasememberAt(co,ip,el) + instance = phase_plasticityInstance(material_phaseAt(co,el)) + call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,instance,of) + case default plasticityType + my_Li = 0.0_pReal + my_dLi_dS = 0.0_pReal + end select plasticityType + + Li = Li + my_Li + dLi_dS = dLi_dS + my_dLi_dS + + KinematicsLoop: do k = 1, phase_Nkinematics(material_phaseAt(co,el)) + kinematicsType: select case (phase_kinematics(k,material_phaseAt(co,el))) + case (KINEMATICS_cleavage_opening_ID) kinematicsType + call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el) + case (KINEMATICS_slipplane_opening_ID) kinematicsType + call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el) + case (KINEMATICS_thermal_expansion_ID) kinematicsType + me = material_phaseMemberAt(co,ip,el) + ph = material_phaseAt(co,el) + call thermalexpansion_LiAndItsTangent(my_Li, my_dLi_dS, ph,me) + case default kinematicsType + my_Li = 0.0_pReal + my_dLi_dS = 0.0_pReal + end select kinematicsType + Li = Li + my_Li + dLi_dS = dLi_dS + my_dLi_dS + enddo KinematicsLoop + + FiInv = math_inv33(Fi) + detFi = math_det33(Fi) + Li = matmul(matmul(Fi,Li),FiInv)*detFi !< push forward to intermediate configuration + temp_33 = matmul(FiInv,Li) + + do i = 1,3; do j = 1,3 + dLi_dS(1:3,1:3,i,j) = matmul(matmul(Fi,dLi_dS(1:3,1:3,i,j)),FiInv)*detFi + dLi_dFi(1:3,1:3,i,j) = dLi_dFi(1:3,1:3,i,j) + Li*FiInv(j,i) + dLi_dFi(1:3,i,1:3,j) = dLi_dFi(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i) + enddo; enddo + +end subroutine constitutive_LiAndItsTangents + + +end submodule eigendeformation diff --git a/src/kinematics_cleavage_opening.f90 b/src/phase_mechanics_eigendeformation_cleavageopening.f90 similarity index 96% rename from src/kinematics_cleavage_opening.f90 rename to src/phase_mechanics_eigendeformation_cleavageopening.f90 index a29a290f8..59be7837a 100644 --- a/src/kinematics_cleavage_opening.f90 +++ b/src/phase_mechanics_eigendeformation_cleavageopening.f90 @@ -4,7 +4,7 @@ !> @brief material subroutine incorporating kinematics resulting from opening of cleavage planes !> @details to be done !-------------------------------------------------------------------------------------------------- -submodule(constitutive:constitutive_damage) kinematics_cleavage_opening +submodule(phase:eigendeformation) cleavageopening integer, dimension(:), allocatable :: kinematics_cleavage_opening_instance @@ -31,8 +31,8 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- module function kinematics_cleavage_opening_init(kinematics_length) result(myKinematics) - - integer, intent(in) :: kinematics_length + + integer, intent(in) :: kinematics_length logical, dimension(:,:), allocatable :: myKinematics integer :: Ninstances,p,k @@ -42,9 +42,9 @@ module function kinematics_cleavage_opening_init(kinematics_length) result(myKin phases, & phase, & kinematics, & - kinematic_type - - print'(/,a)', ' <<<+- kinematics_cleavage_opening init -+>>>' + kinematic_type + + print'(/,a)', ' <<<+- phase:mechanics:eigendeformation:cleavageopening init -+>>>' myKinematics = kinematics_active('cleavage_opening',kinematics_length) Ninstances = count(myKinematics) @@ -63,7 +63,7 @@ module function kinematics_cleavage_opening_init(kinematics_length) result(myKin do k = 1, kinematics%length if(myKinematics(k,p)) then associate(prm => param(kinematics_cleavage_opening_instance(p))) - kinematic_type => kinematics%get(k) + kinematic_type => kinematics%get(k) N_cl = kinematic_type%get_asInts('N_cl') prm%sum_N_cl = sum(abs(N_cl)) @@ -162,4 +162,4 @@ module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, end subroutine kinematics_cleavage_opening_LiAndItsTangent -end submodule kinematics_cleavage_opening +end submodule cleavageopening diff --git a/src/kinematics_slipplane_opening.f90 b/src/phase_mechanics_eigendeformation_slipplaneopening.f90 similarity index 96% rename from src/kinematics_slipplane_opening.f90 rename to src/phase_mechanics_eigendeformation_slipplaneopening.f90 index 84edab122..a144d39a6 100644 --- a/src/kinematics_slipplane_opening.f90 +++ b/src/phase_mechanics_eigendeformation_slipplaneopening.f90 @@ -4,7 +4,7 @@ !> @brief material subroutine incorporating kinematics resulting from opening of slip planes !> @details to be done !-------------------------------------------------------------------------------------------------- -submodule(constitutive:constitutive_damage) kinematics_slipplane_opening +submodule(phase:eigendeformation) slipplaneopening integer, dimension(:), allocatable :: kinematics_slipplane_opening_instance @@ -34,7 +34,7 @@ contains !-------------------------------------------------------------------------------------------------- module function kinematics_slipplane_opening_init(kinematics_length) result(myKinematics) - integer, intent(in) :: kinematics_length + integer, intent(in) :: kinematics_length logical, dimension(:,:), allocatable :: myKinematics integer :: Ninstances,p,i,k @@ -47,9 +47,9 @@ module function kinematics_slipplane_opening_init(kinematics_length) result(myKi mech, & pl, & kinematics, & - kinematic_type - - print'(/,a)', ' <<<+- kinematics_slipplane init -+>>>' + kinematic_type + + print'(/,a)', ' <<<+- phase:mechanics:eigendeformation:slipplaneopening init -+>>>' myKinematics = kinematics_active('slipplane_opening',kinematics_length) Ninstances = count(myKinematics) @@ -70,7 +70,7 @@ module function kinematics_slipplane_opening_init(kinematics_length) result(myKi do k = 1, kinematics%length if(myKinematics(k,p)) then associate(prm => param(kinematics_slipplane_opening_instance(p))) - kinematic_type => kinematics%get(k) + kinematic_type => kinematics%get(k) prm%dot_o = kinematic_type%get_asFloat('dot_o') prm%q = kinematic_type%get_asFloat('q') @@ -193,4 +193,4 @@ module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S end subroutine kinematics_slipplane_opening_LiAndItsTangent -end submodule kinematics_slipplane_opening +end submodule slipplaneopening diff --git a/src/kinematics_thermal_expansion.f90 b/src/phase_mechanics_eigendeformation_thermalexpansion.f90 similarity index 78% rename from src/kinematics_thermal_expansion.f90 rename to src/phase_mechanics_eigendeformation_thermalexpansion.f90 index 6e4fa8263..6630f0eb7 100644 --- a/src/kinematics_thermal_expansion.f90 +++ b/src/phase_mechanics_eigendeformation_thermalexpansion.f90 @@ -3,7 +3,7 @@ !> @brief material subroutine incorporating kinematics resulting from thermal expansion !> @details to be done !-------------------------------------------------------------------------------------------------- -submodule(constitutive:constitutive_thermal) kinematics_thermal_expansion +submodule(phase:eigendeformation) thermalexpansion integer, dimension(:), allocatable :: kinematics_thermal_expansion_instance @@ -16,7 +16,6 @@ submodule(constitutive:constitutive_thermal) kinematics_thermal_expansion type(tParameters), dimension(:), allocatable :: param - contains @@ -37,7 +36,7 @@ module function kinematics_thermal_expansion_init(kinematics_length) result(myKi kinematics, & kinematic_type - print'(/,a)', ' <<<+- kinematics_thermal_expansion init -+>>>' + print'(/,a)', ' <<<+- phase:mechanics:eigendeformation:thermalexpansion init -+>>>' myKinematics = kinematics_active('thermal_expansion',kinematics_length) Ninstances = count(myKinematics) @@ -84,7 +83,7 @@ end function kinematics_thermal_expansion_init !-------------------------------------------------------------------------------------------------- !> @brief constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ph,me) +module subroutine thermalexpansion_LiAndItsTangent(Li, dLi_dTstar, ph,me) integer, intent(in) :: ph, me real(pReal), intent(out), dimension(3,3) :: & @@ -92,28 +91,25 @@ module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, p real(pReal), intent(out), dimension(3,3,3,3) :: & dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero) - integer :: & - phase, & - homog real(pReal) :: T, dot_T - T = current(ph)%T(me) - dot_T = current(ph)%dot_T(me) + T = thermal_T(ph,me) + dot_T = thermal_dot_T(ph,me) associate(prm => param(kinematics_thermal_expansion_instance(ph))) - Li = dot_T * ( & - prm%A(1:3,1:3,1)*(T - prm%T_ref)**0 & ! constant coefficient - + prm%A(1:3,1:3,2)*(T - prm%T_ref)**1 & ! linear coefficient - + prm%A(1:3,1:3,3)*(T - prm%T_ref)**2 & ! quadratic coefficient - ) / & - (1.0_pReal & - + prm%A(1:3,1:3,1)*(T - prm%T_ref)**1 / 1. & - + prm%A(1:3,1:3,2)*(T - prm%T_ref)**2 / 2. & - + prm%A(1:3,1:3,3)*(T - prm%T_ref)**3 / 3. & - ) - end associate + Li = dot_T * ( & + prm%A(1:3,1:3,1)*(T - prm%T_ref)**0 & ! constant coefficient + + prm%A(1:3,1:3,2)*(T - prm%T_ref)**1 & ! linear coefficient + + prm%A(1:3,1:3,3)*(T - prm%T_ref)**2 & ! quadratic coefficient + ) / & + (1.0_pReal & + + prm%A(1:3,1:3,1)*(T - prm%T_ref)**1 / 1. & + + prm%A(1:3,1:3,2)*(T - prm%T_ref)**2 / 2. & + + prm%A(1:3,1:3,3)*(T - prm%T_ref)**3 / 3. & + ) + end associate dLi_dTstar = 0.0_pReal -end subroutine kinematics_thermal_expansion_LiAndItsTangent +end subroutine thermalexpansion_LiAndItsTangent -end submodule kinematics_thermal_expansion +end submodule thermalexpansion diff --git a/src/phase_mechanics_plastic.f90 b/src/phase_mechanics_plastic.f90 new file mode 100644 index 000000000..fff1fbdcf --- /dev/null +++ b/src/phase_mechanics_plastic.f90 @@ -0,0 +1,472 @@ +submodule(phase:mechanics) plastic + + interface + + module function plastic_none_init() result(myPlasticity) + logical, dimension(:), allocatable :: & + myPlasticity + end function plastic_none_init + + module function plastic_isotropic_init() result(myPlasticity) + logical, dimension(:), allocatable :: & + myPlasticity + end function plastic_isotropic_init + + module function plastic_phenopowerlaw_init() result(myPlasticity) + logical, dimension(:), allocatable :: & + myPlasticity + end function plastic_phenopowerlaw_init + + module function plastic_kinehardening_init() result(myPlasticity) + logical, dimension(:), allocatable :: & + myPlasticity + end function plastic_kinehardening_init + + module function plastic_dislotwin_init() result(myPlasticity) + logical, dimension(:), allocatable :: & + myPlasticity + end function plastic_dislotwin_init + + module function plastic_dislotungsten_init() result(myPlasticity) + logical, dimension(:), allocatable :: & + myPlasticity + end function plastic_dislotungsten_init + + module function plastic_nonlocal_init() result(myPlasticity) + logical, dimension(:), allocatable :: & + myPlasticity + end function plastic_nonlocal_init + + module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) + real(pReal), dimension(3,3), intent(out) :: & + Lp + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp + real(pReal), dimension(3,3), intent(in) :: & + Mp + integer, intent(in) :: & + ph, & + me + end subroutine isotropic_LpAndItsTangent + + pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) + real(pReal), dimension(3,3), intent(out) :: & + Lp + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp + real(pReal), dimension(3,3), intent(in) :: & + Mp + integer, intent(in) :: & + ph, & + me + end subroutine phenopowerlaw_LpAndItsTangent + + pure module subroutine kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) + real(pReal), dimension(3,3), intent(out) :: & + Lp + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp + real(pReal), dimension(3,3), intent(in) :: & + Mp + integer, intent(in) :: & + ph, & + me + end subroutine kinehardening_LpAndItsTangent + + module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me) + real(pReal), dimension(3,3), intent(out) :: & + Lp + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp + + real(pReal), dimension(3,3), intent(in) :: & + Mp + real(pReal), intent(in) :: & + T + integer, intent(in) :: & + ph, & + me + end subroutine dislotwin_LpAndItsTangent + + pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me) + real(pReal), dimension(3,3), intent(out) :: & + Lp + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp + + real(pReal), dimension(3,3), intent(in) :: & + Mp + real(pReal), intent(in) :: & + T + integer, intent(in) :: & + ph, & + me + end subroutine dislotungsten_LpAndItsTangent + + module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & + Mp,Temperature,ph,me,ip,el) + real(pReal), dimension(3,3), intent(out) :: & + Lp + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + Temperature + integer, intent(in) :: & + ph, & + me, & + ip, & !< current integration point + el !< current element number + end subroutine nonlocal_LpAndItsTangent + + + module subroutine isotropic_dotState(Mp,ph,me) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + ph, & + me + end subroutine isotropic_dotState + + module subroutine phenopowerlaw_dotState(Mp,ph,me) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + ph, & + me + end subroutine phenopowerlaw_dotState + + module subroutine plastic_kinehardening_dotState(Mp,ph,me) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + ph, & + me + end subroutine plastic_kinehardening_dotState + + module subroutine dislotwin_dotState(Mp,T,ph,me) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + T + integer, intent(in) :: & + ph, & + me + end subroutine dislotwin_dotState + + module subroutine dislotungsten_dotState(Mp,T,ph,me) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + T + integer, intent(in) :: & + ph, & + me + end subroutine dislotungsten_dotState + + module subroutine nonlocal_dotState(Mp,Temperature,timestep,ph,me,ip,el) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< MandelStress + real(pReal), intent(in) :: & + Temperature, & !< temperature + timestep !< substepped crystallite time increment + integer, intent(in) :: & + ph, & + me, & + ip, & !< current integration point + el !< current element number + end subroutine nonlocal_dotState + + module subroutine dislotwin_dependentState(T,instance,me) + integer, intent(in) :: & + instance, & + me + real(pReal), intent(in) :: & + T + end subroutine dislotwin_dependentState + + module subroutine dislotungsten_dependentState(instance,me) + integer, intent(in) :: & + instance, & + me + end subroutine dislotungsten_dependentState + + module subroutine nonlocal_dependentState(instance, me, ip, el) + integer, intent(in) :: & + instance, & + me, & + ip, & !< current integration point + el !< current element number + end subroutine nonlocal_dependentState + + module subroutine plastic_kinehardening_deltaState(Mp,instance,me) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + me + end subroutine plastic_kinehardening_deltaState + + module subroutine plastic_nonlocal_deltaState(Mp,instance,me,ip,el) + real(pReal), dimension(3,3), intent(in) :: & + Mp + integer, intent(in) :: & + instance, & + me, & + ip, & + el + end subroutine plastic_nonlocal_deltaState + + end interface + +contains + +module subroutine plastic_init + + + print'(/,a)', ' <<<+- phase:mechanics:plastic init -+>>>' + + where(plastic_none_init()) phase_plasticity = PLASTICITY_NONE_ID + where(plastic_isotropic_init()) phase_plasticity = PLASTICITY_ISOTROPIC_ID + where(plastic_phenopowerlaw_init()) phase_plasticity = PLASTICITY_PHENOPOWERLAW_ID + where(plastic_kinehardening_init()) phase_plasticity = PLASTICITY_KINEHARDENING_ID + where(plastic_dislotwin_init()) phase_plasticity = PLASTICITY_DISLOTWIN_ID + where(plastic_dislotungsten_init()) phase_plasticity = PLASTICITY_DISLOTUNGSTEN_ID + where(plastic_nonlocal_init()) phase_plasticity = PLASTICITY_NONLOCAL_ID + +end subroutine plastic_init + +!-------------------------------------------------------------------------------------------------- +!> @brief contains the constitutive equation for calculating the velocity gradient +! ToDo: Discuss whether it makes sense if crystallite handles the configuration conversion, i.e. +! Mp in, dLp_dMp out +!-------------------------------------------------------------------------------------------------- +module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & + S, Fi, co, ip, el) + integer, intent(in) :: & + co, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + S, & !< 2nd Piola-Kirchhoff stress + Fi !< intermediate deformation gradient + real(pReal), intent(out), dimension(3,3) :: & + Lp !< plastic velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLp_dS, & + dLp_dFi !< derivative me Lp with respect to Fi + + real(pReal), dimension(3,3,3,3) :: & + dLp_dMp !< derivative of Lp with respect to Mandel stress + real(pReal), dimension(3,3) :: & + Mp !< Mandel stress work conjugate with Lp + integer :: & + i, j, me, ph + + + Mp = matmul(matmul(transpose(Fi),Fi),S) + me = material_phasememberAt(co,ip,el) + ph = material_phaseAt(co,el) + + plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) + + case (PLASTICITY_NONE_ID) plasticityType + Lp = 0.0_pReal + dLp_dMp = 0.0_pReal + + case (PLASTICITY_ISOTROPIC_ID) plasticityType + call isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) + + case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType + call phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) + + case (PLASTICITY_KINEHARDENING_ID) plasticityType + call kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) + + case (PLASTICITY_NONLOCAL_ID) plasticityType + call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me,ip,el) + + case (PLASTICITY_DISLOTWIN_ID) plasticityType + call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me) + + case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType + call dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me) + + end select plasticityType + + do i=1,3; do j=1,3 + dLp_dFi(i,j,1:3,1:3) = matmul(matmul(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + & + matmul(matmul(Fi,dLp_dMp(i,j,1:3,1:3)),S) + dLp_dS(i,j,1:3,1:3) = matmul(matmul(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi) + enddo; enddo + +end subroutine plastic_LpAndItsTangents + + +!-------------------------------------------------------------------------------------------------- +!> @brief contains the constitutive equation for calculating the rate of change of microstructure +!-------------------------------------------------------------------------------------------------- +module function plastic_dotState(subdt,co,ip,el,ph,me) result(broken) + + integer, intent(in) :: & + co, & !< component-ID of integration point + ip, & !< integration point + el, & !< element + ph, & + me + real(pReal), intent(in) :: & + subdt !< timestep + real(pReal), dimension(3,3) :: & + Mp + logical :: broken + + + Mp = matmul(matmul(transpose(constitutive_mech_Fi(ph)%data(1:3,1:3,me)),& + constitutive_mech_Fi(ph)%data(1:3,1:3,me)),constitutive_mech_S(ph)%data(1:3,1:3,me)) + + plasticityType: select case (phase_plasticity(ph)) + + case (PLASTICITY_ISOTROPIC_ID) plasticityType + call isotropic_dotState(Mp,ph,me) + + case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType + call phenopowerlaw_dotState(Mp,ph,me) + + case (PLASTICITY_KINEHARDENING_ID) plasticityType + call plastic_kinehardening_dotState(Mp,ph,me) + + case (PLASTICITY_DISLOTWIN_ID) plasticityType + call dislotwin_dotState(Mp,thermal_T(ph,me),ph,me) + + case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType + call dislotungsten_dotState(Mp,thermal_T(ph,me),ph,me) + + case (PLASTICITY_NONLOCAL_ID) plasticityType + call nonlocal_dotState(Mp,thermal_T(ph,me),subdt,ph,me,ip,el) + end select plasticityType + broken = any(IEEE_is_NaN(plasticState(ph)%dotState(:,me))) + + +end function plastic_dotState + + +!-------------------------------------------------------------------------------------------------- +!> @brief calls microstructure function of the different plasticity constitutive models +!-------------------------------------------------------------------------------------------------- +module subroutine plastic_dependentState(co, ip, el) + + integer, intent(in) :: & + co, & !< component-ID of integration point + ip, & !< integration point + el !< element + + integer :: & + ph, & + instance, me + + + ph = material_phaseAt(co,el) + me = material_phasememberAt(co,ip,el) + instance = phase_plasticityInstance(ph) + + plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) + + case (PLASTICITY_DISLOTWIN_ID) plasticityType + call dislotwin_dependentState(thermal_T(ph,me),instance,me) + + case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType + call dislotungsten_dependentState(instance,me) + + case (PLASTICITY_NONLOCAL_ID) plasticityType + call nonlocal_dependentState(instance,me,ip,el) + + end select plasticityType + +end subroutine plastic_dependentState + + +!-------------------------------------------------------------------------------------------------- +!> @brief for constitutive models having an instantaneous change of state +!> will return false if delta state is not needed/supported by the constitutive model +!-------------------------------------------------------------------------------------------------- +module function plastic_deltaState(co, ip, el, ph, me) result(broken) + + integer, intent(in) :: & + co, & !< component-ID of integration point + ip, & !< integration point + el, & !< element + ph, & + me + logical :: & + broken + + real(pReal), dimension(3,3) :: & + Mp + integer :: & + instance, & + myOffset, & + mySize + + + Mp = matmul(matmul(transpose(constitutive_mech_Fi(ph)%data(1:3,1:3,me)),& + constitutive_mech_Fi(ph)%data(1:3,1:3,me)),constitutive_mech_S(ph)%data(1:3,1:3,me)) + instance = phase_plasticityInstance(ph) + + plasticityType: select case (phase_plasticity(ph)) + + case (PLASTICITY_KINEHARDENING_ID) plasticityType + call plastic_kinehardening_deltaState(Mp,instance,me) + broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,me))) + + case (PLASTICITY_NONLOCAL_ID) plasticityType + call plastic_nonlocal_deltaState(Mp,instance,me,ip,el) + broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,me))) + + case default + broken = .false. + + end select plasticityType + + if(.not. broken) then + select case(phase_plasticity(ph)) + case (PLASTICITY_NONLOCAL_ID,PLASTICITY_KINEHARDENING_ID) + + myOffset = plasticState(ph)%offsetDeltaState + mySize = plasticState(ph)%sizeDeltaState + plasticState(ph)%state(myOffset + 1:myOffset + mySize,me) = & + plasticState(ph)%state(myOffset + 1:myOffset + mySize,me) + plasticState(ph)%deltaState(1:mySize,me) + end select + endif + +end function plastic_deltaState + + +!-------------------------------------------------------------------------------------------------- +!> @brief checks if a plastic module is active or not +!-------------------------------------------------------------------------------------------------- +function plastic_active(plastic_label) result(active_plastic) + + character(len=*), intent(in) :: plastic_label !< type of plasticity model + logical, dimension(:), allocatable :: active_plastic + + class(tNode), pointer :: & + phases, & + phase, & + mech, & + pl + integer :: ph + + phases => config_material%get('phase') + allocate(active_plastic(phases%length), source = .false. ) + do ph = 1, phases%length + phase => phases%get(ph) + mech => phase%get('mechanics') + pl => mech%get('plasticity') + if(pl%get_asString('type') == plastic_label) active_plastic(ph) = .true. + enddo + +end function plastic_active + +end submodule plastic diff --git a/src/constitutive_plastic_disloTungsten.f90 b/src/phase_mechanics_plastic_dislotungsten.f90 similarity index 91% rename from src/constitutive_plastic_disloTungsten.f90 rename to src/phase_mechanics_plastic_dislotungsten.f90 index c39ae5c2b..a47132c63 100644 --- a/src/constitutive_plastic_disloTungsten.f90 +++ b/src/phase_mechanics_plastic_dislotungsten.f90 @@ -5,7 +5,7 @@ !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief crystal plasticity model for bcc metals, especially Tungsten !-------------------------------------------------------------------------------------------------- -submodule(constitutive:constitutive_mech) plastic_dislotungsten +submodule(phase:plastic) dislotungsten real(pReal), parameter :: & kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin @@ -17,7 +17,7 @@ submodule(constitutive:constitutive_mech) plastic_dislotungsten D_0 = 1.0_pReal, & !< prefactor for self-diffusion coefficient Q_cl = 1.0_pReal !< activation energy for dislocation climb real(pReal), allocatable, dimension(:) :: & - b_sl, & !< magnitude of Burgers vector [m] + b_sl, & !< magnitude me Burgers vector [m] D_a, & i_sl, & !< Adj. parameter for distance between 2 forest dislocations f_at, & !< factor to calculate atomic volume @@ -97,13 +97,13 @@ module function plastic_dislotungsten_init() result(myPlasticity) mech, & pl - print'(/,a)', ' <<<+- plastic_dislotungsten init -+>>>' + print'(/,a)', ' <<<+- phase:mechanics:plastic:dislotungsten init -+>>>' myPlasticity = plastic_active('dislotungsten') Ninstances = count(myPlasticity) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) if(Ninstances == 0) return - + print*, 'Cereceda et al., International Journal of Plasticity 78:242–256, 2016' print*, 'https://dx.doi.org/10.1016/j.ijplas.2015.09.002' @@ -222,7 +222,7 @@ module function plastic_dislotungsten_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! allocate state arrays - Nconstituents = count(material_phaseAt == p) * discretization_nIPs + Nconstituents = count(material_phaseAt2 == p) sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl sizeState = sizeDotState @@ -272,8 +272,8 @@ end function plastic_dislotungsten_init !-------------------------------------------------------------------------------------------------- !> @brief Calculate plastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- -pure module subroutine plastic_dislotungsten_LpAndItsTangent(Lp,dLp_dMp, & - Mp,T,instance,of) +pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp, & + Mp,T,ph,me) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & @@ -284,21 +284,21 @@ pure module subroutine plastic_dislotungsten_LpAndItsTangent(Lp,dLp_dMp, & real(pReal), intent(in) :: & T !< temperature integer, intent(in) :: & - instance, & - of + ph, & + me integer :: & i,k,l,m,n - real(pReal), dimension(param(instance)%sum_N_sl) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & dot_gamma_pos,dot_gamma_neg, & ddot_gamma_dtau_pos,ddot_gamma_dtau_neg Lp = 0.0_pReal dLp_dMp = 0.0_pReal - associate(prm => param(instance)) + associate(prm => param(phase_plasticityInstance(ph))) - call kinetics(Mp,T,instance,of,dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg) + call kinetics(Mp,T,phase_plasticityInstance(ph),me,dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg) do i = 1, prm%sum_N_sl Lp = Lp + (dot_gamma_pos(i)+dot_gamma_neg(i))*prm%P_sl(1:3,1:3,i) forall (k=1:3,l=1:3,m=1:3,n=1:3) & @@ -309,25 +309,25 @@ pure module subroutine plastic_dislotungsten_LpAndItsTangent(Lp,dLp_dMp, & end associate -end subroutine plastic_dislotungsten_LpAndItsTangent +end subroutine dislotungsten_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- -module subroutine plastic_dislotungsten_dotState(Mp,T,instance,of) +module subroutine dislotungsten_dotState(Mp,T,ph,me) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress real(pReal), intent(in) :: & T !< temperature integer, intent(in) :: & - instance, & - of + ph, & + me real(pReal) :: & VacancyDiffusion - real(pReal), dimension(param(instance)%sum_N_sl) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & gdot_pos, gdot_neg,& tau_pos,& tau_neg, & @@ -336,13 +336,14 @@ module subroutine plastic_dislotungsten_dotState(Mp,T,instance,of) dot_rho_dip_climb, & dip_distance - associate(prm => param(instance), stt => state(instance),dot => dotState(instance), dst => dependentState(instance)) + associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph)),& + dot => dotState(phase_plasticityInstance(ph)), dst => dependentState(phase_plasticityInstance(ph))) - call kinetics(Mp,T,instance,of,& + call kinetics(Mp,T,phase_plasticityInstance(ph),me,& gdot_pos,gdot_neg, & tau_pos_out = tau_pos,tau_neg_out = tau_neg) - dot%gamma_sl(:,of) = (gdot_pos+gdot_neg) ! ToDo: needs to be abs + dot%gamma_sl(:,me) = (gdot_pos+gdot_neg) ! ToDo: needs to be abs VacancyDiffusion = prm%D_0*exp(-prm%Q_cl/(kB*T)) where(dEq0(tau_pos)) ! ToDo: use avg of pos and neg @@ -351,50 +352,50 @@ module subroutine plastic_dislotungsten_dotState(Mp,T,instance,of) else where dip_distance = math_clip(3.0_pReal*prm%mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos)), & prm%D_a, & ! lower limit - dst%Lambda_sl(:,of)) ! upper limit - dot_rho_dip_formation = merge(2.0_pReal*dip_distance* stt%rho_mob(:,of)*abs(dot%gamma_sl(:,of))/prm%b_sl, & ! ToDo: ignore region of spontaneous annihilation + dst%Lambda_sl(:,me)) ! upper limit + dot_rho_dip_formation = merge(2.0_pReal*dip_distance* stt%rho_mob(:,me)*abs(dot%gamma_sl(:,me))/prm%b_sl, & ! ToDo: ignore region of spontaneous annihilation 0.0_pReal, & prm%dipoleformation) v_cl = (3.0_pReal*prm%mu*VacancyDiffusion*prm%f_at/(2.0_pReal*pi*kB*T)) & * (1.0_pReal/(dip_distance+prm%D_a)) - dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,of))/(dip_distance-prm%D_a) ! ToDo: Discuss with Franz: Stress dependency? + dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,me))/(dip_distance-prm%D_a) ! ToDo: Discuss with Franz: Stress dependency? end where - dot%rho_mob(:,of) = abs(dot%gamma_sl(:,of))/(prm%b_sl*dst%Lambda_sl(:,of)) & ! multiplication + dot%rho_mob(:,me) = abs(dot%gamma_sl(:,me))/(prm%b_sl*dst%Lambda_sl(:,me)) & ! multiplication - dot_rho_dip_formation & - - (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_mob(:,of)*abs(dot%gamma_sl(:,of)) ! Spontaneous annihilation of 2 single edge dislocations - dot%rho_dip(:,of) = dot_rho_dip_formation & - - (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_dip(:,of)*abs(dot%gamma_sl(:,of)) & ! Spontaneous annihilation of a single edge dislocation with a dipole constituent + - (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_mob(:,me)*abs(dot%gamma_sl(:,me)) ! Spontaneous annihilation of 2 single edge dislocations + dot%rho_dip(:,me) = dot_rho_dip_formation & + - (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_dip(:,me)*abs(dot%gamma_sl(:,me)) & ! Spontaneous annihilation of a single edge dislocation with a dipole constituent - dot_rho_dip_climb end associate -end subroutine plastic_dislotungsten_dotState +end subroutine dislotungsten_dotState !-------------------------------------------------------------------------------------------------- !> @brief Calculate derived quantities from state. !-------------------------------------------------------------------------------------------------- -module subroutine plastic_dislotungsten_dependentState(instance,of) +module subroutine dislotungsten_dependentState(instance,me) integer, intent(in) :: & instance, & - of + me real(pReal), dimension(param(instance)%sum_N_sl) :: & dislocationSpacing associate(prm => param(instance), stt => state(instance),dst => dependentState(instance)) - dislocationSpacing = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,of)+stt%rho_dip(:,of))) - dst%threshold_stress(:,of) = prm%mu*prm%b_sl & - * sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,of)+stt%rho_dip(:,of))) + dislocationSpacing = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,me)+stt%rho_dip(:,me))) + dst%threshold_stress(:,me) = prm%mu*prm%b_sl & + * sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,me)+stt%rho_dip(:,me))) - dst%Lambda_sl(:,of) = prm%D/(1.0_pReal+prm%D*dislocationSpacing/prm%i_sl) + dst%Lambda_sl(:,me) = prm%D/(1.0_pReal+prm%D*dislocationSpacing/prm%i_sl) end associate -end subroutine plastic_dislotungsten_dependentState +end subroutine dislotungsten_dependentState !-------------------------------------------------------------------------------------------------- @@ -439,7 +440,7 @@ end subroutine plastic_dislotungsten_results ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end !-------------------------------------------------------------------------------------------------- -pure subroutine kinetics(Mp,T,instance,of, & +pure subroutine kinetics(Mp,T,instance,me, & dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg,tau_pos_out,tau_neg_out) real(pReal), dimension(3,3), intent(in) :: & @@ -448,7 +449,7 @@ pure subroutine kinetics(Mp,T,instance,of, & T !< temperature integer, intent(in) :: & instance, & - of + me real(pReal), intent(out), dimension(param(instance)%sum_N_sl) :: & dot_gamma_pos, & @@ -479,11 +480,11 @@ pure subroutine kinetics(Mp,T,instance,of, & if (present(tau_neg_out)) tau_neg_out = tau_neg associate(BoltzmannRatio => prm%Q_s/(kB*T), & - dot_gamma_0 => stt%rho_mob(:,of)*prm%b_sl*prm%v_0, & - effectiveLength => dst%Lambda_sl(:,of) - prm%w) + dot_gamma_0 => stt%rho_mob(:,me)*prm%b_sl*prm%v_0, & + effectiveLength => dst%Lambda_sl(:,me) - prm%w) - significantPositiveTau: where(abs(tau_pos)-dst%threshold_stress(:,of) > tol_math_check) - StressRatio = (abs(tau_pos)-dst%threshold_stress(:,of))/prm%tau_Peierls + significantPositiveTau: where(abs(tau_pos)-dst%threshold_stress(:,me) > tol_math_check) + StressRatio = (abs(tau_pos)-dst%threshold_stress(:,me))/prm%tau_Peierls StressRatio_p = StressRatio** prm%p StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) @@ -499,7 +500,7 @@ pure subroutine kinetics(Mp,T,instance,of, & end where significantPositiveTau if (present(ddot_gamma_dtau_pos)) then - significantPositiveTau2: where(abs(tau_pos)-dst%threshold_stress(:,of) > tol_math_check) + significantPositiveTau2: where(abs(tau_pos)-dst%threshold_stress(:,me) > tol_math_check) dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_Peierls dtk = -1.0_pReal * t_k / tau_pos @@ -512,8 +513,8 @@ pure subroutine kinetics(Mp,T,instance,of, & end where significantPositiveTau2 endif - significantNegativeTau: where(abs(tau_neg)-dst%threshold_stress(:,of) > tol_math_check) - StressRatio = (abs(tau_neg)-dst%threshold_stress(:,of))/prm%tau_Peierls + significantNegativeTau: where(abs(tau_neg)-dst%threshold_stress(:,me) > tol_math_check) + StressRatio = (abs(tau_neg)-dst%threshold_stress(:,me))/prm%tau_Peierls StressRatio_p = StressRatio** prm%p StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) @@ -529,7 +530,7 @@ pure subroutine kinetics(Mp,T,instance,of, & end where significantNegativeTau if (present(ddot_gamma_dtau_neg)) then - significantNegativeTau2: where(abs(tau_neg)-dst%threshold_stress(:,of) > tol_math_check) + significantNegativeTau2: where(abs(tau_neg)-dst%threshold_stress(:,me) > tol_math_check) dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_Peierls dtk = -1.0_pReal * t_k / tau_neg @@ -547,4 +548,4 @@ pure subroutine kinetics(Mp,T,instance,of, & end subroutine kinetics -end submodule plastic_dislotungsten +end submodule dislotungsten diff --git a/src/constitutive_plastic_dislotwin.f90 b/src/phase_mechanics_plastic_dislotwin.f90 similarity index 91% rename from src/constitutive_plastic_dislotwin.f90 rename to src/phase_mechanics_plastic_dislotwin.f90 index dea84f751..9f7464323 100644 --- a/src/constitutive_plastic_dislotwin.f90 +++ b/src/phase_mechanics_plastic_dislotwin.f90 @@ -7,7 +7,7 @@ !> @brief material subroutine incoprorating dislocation and twinning physics !> @details to be done !-------------------------------------------------------------------------------------------------- -submodule(constitutive:constitutive_mech) plastic_dislotwin +submodule(phase:plastic) dislotwin real(pReal), parameter :: & kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin @@ -144,7 +144,7 @@ module function plastic_dislotwin_init() result(myPlasticity) mech, & pl - print'(/,a)', ' <<<+- plastic_dislotwin init -+>>>' + print'(/,a)', ' <<<+- phase:mechanics:plastic:dislotwin init -+>>>' myPlasticity = plastic_active('dislotwin') Ninstances = count(myPlasticity) @@ -408,7 +408,7 @@ module function plastic_dislotwin_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! allocate state arrays - Nconstituents = count(material_phaseAt == p) * discretization_nIPs + Nconstituents = count(material_phaseAt2 == p) sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl & + size(['f_tw']) * prm%sum_N_tw & + size(['f_tr']) * prm%sum_N_tr @@ -521,12 +521,12 @@ end function plastic_dislotwin_homogenizedC !-------------------------------------------------------------------------------------------------- !> @brief Calculate plastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- -module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) +module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me) real(pReal), dimension(3,3), intent(out) :: Lp real(pReal), dimension(3,3,3,3), intent(out) :: dLp_dMp real(pReal), dimension(3,3), intent(in) :: Mp - integer, intent(in) :: instance,of + integer, intent(in) :: ph,me real(pReal), intent(in) :: T integer :: i,k,l,m,n @@ -535,11 +535,11 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) BoltzmannRatio, & ddot_gamma_dtau, & tau - real(pReal), dimension(param(instance)%sum_N_sl) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & dot_gamma_sl,ddot_gamma_dtau_slip - real(pReal), dimension(param(instance)%sum_N_tw) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_tw) :: & dot_gamma_twin,ddot_gamma_dtau_twin - real(pReal), dimension(param(instance)%sum_N_tr) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_tr) :: & dot_gamma_tr,ddot_gamma_dtau_trans real(pReal):: dot_gamma_sb real(pReal), dimension(3,3) :: eigVectors, P_sb @@ -564,16 +564,16 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) 0, 1, 1 & ],pReal),[ 3,6]) - associate(prm => param(instance), stt => state(instance)) + associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph))) f_unrotated = 1.0_pReal & - - sum(stt%f_tw(1:prm%sum_N_tw,of)) & - - sum(stt%f_tr(1:prm%sum_N_tr,of)) + - sum(stt%f_tw(1:prm%sum_N_tw,me)) & + - sum(stt%f_tr(1:prm%sum_N_tr,me)) Lp = 0.0_pReal dLp_dMp = 0.0_pReal - call kinetics_slip(Mp,T,instance,of,dot_gamma_sl,ddot_gamma_dtau_slip) + call kinetics_slip(Mp,T,phase_plasticityInstance(ph),me,dot_gamma_sl,ddot_gamma_dtau_slip) slipContribution: do i = 1, prm%sum_N_sl Lp = Lp + dot_gamma_sl(i)*prm%P_sl(1:3,1:3,i) forall (k=1:3,l=1:3,m=1:3,n=1:3) & @@ -581,7 +581,7 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) + ddot_gamma_dtau_slip(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i) enddo slipContribution - call kinetics_twin(Mp,T,dot_gamma_sl,instance,of,dot_gamma_twin,ddot_gamma_dtau_twin) + call kinetics_twin(Mp,T,dot_gamma_sl,phase_plasticityInstance(ph),me,dot_gamma_twin,ddot_gamma_dtau_twin) twinContibution: do i = 1, prm%sum_N_tw Lp = Lp + dot_gamma_twin(i)*prm%P_tw(1:3,1:3,i) forall (k=1:3,l=1:3,m=1:3,n=1:3) & @@ -589,7 +589,7 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) + ddot_gamma_dtau_twin(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i) enddo twinContibution - call kinetics_trans(Mp,T,dot_gamma_sl,instance,of,dot_gamma_tr,ddot_gamma_dtau_trans) + call kinetics_trans(Mp,T,dot_gamma_sl,phase_plasticityInstance(ph),me,dot_gamma_tr,ddot_gamma_dtau_trans) transContibution: do i = 1, prm%sum_N_tr Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i) forall (k=1:3,l=1:3,m=1:3,n=1:3) & @@ -628,21 +628,21 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) end associate -end subroutine plastic_dislotwin_LpAndItsTangent +end subroutine dislotwin_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- -module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) +module subroutine dislotwin_dotState(Mp,T,ph,me) real(pReal), dimension(3,3), intent(in):: & Mp !< Mandel stress real(pReal), intent(in) :: & T !< temperature at integration point integer, intent(in) :: & - instance, & - of + ph, & + me integer :: i real(pReal) :: & @@ -653,25 +653,25 @@ module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) tau, & sigma_cl, & !< climb stress b_d !< ratio of Burgers vector to stacking fault width - real(pReal), dimension(param(instance)%sum_N_sl) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & dot_rho_dip_formation, & dot_rho_dip_climb, & rho_dip_distance_min, & dot_gamma_sl - real(pReal), dimension(param(instance)%sum_N_tw) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_tw) :: & dot_gamma_twin - real(pReal), dimension(param(instance)%sum_N_tr) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_tr) :: & dot_gamma_tr - associate(prm => param(instance), stt => state(instance), & - dot => dotState(instance), dst => dependentState(instance)) + associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph)), & + dot => dotState(phase_plasticityInstance(ph)), dst => dependentState(phase_plasticityInstance(ph))) f_unrotated = 1.0_pReal & - - sum(stt%f_tw(1:prm%sum_N_tw,of)) & - - sum(stt%f_tr(1:prm%sum_N_tr,of)) + - sum(stt%f_tw(1:prm%sum_N_tw,me)) & + - sum(stt%f_tr(1:prm%sum_N_tr,me)) - call kinetics_slip(Mp,T,instance,of,dot_gamma_sl) - dot%gamma_sl(:,of) = abs(dot_gamma_sl) + call kinetics_slip(Mp,T,phase_plasticityInstance(ph),me,dot_gamma_sl) + dot%gamma_sl(:,me) = abs(dot_gamma_sl) rho_dip_distance_min = prm%D_a*prm%b_sl @@ -683,12 +683,12 @@ module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) dot_rho_dip_climb(i) = 0.0_pReal else significantSlipStress rho_dip_distance = 3.0_pReal*prm%mu*prm%b_sl(i)/(16.0_pReal*PI*abs(tau)) - rho_dip_distance = math_clip(rho_dip_distance, right = dst%Lambda_sl(i,of)) + rho_dip_distance = math_clip(rho_dip_distance, right = dst%Lambda_sl(i,me)) rho_dip_distance = math_clip(rho_dip_distance, left = rho_dip_distance_min(i)) if (prm%dipoleFormation) then dot_rho_dip_formation(i) = 2.0_pReal*(rho_dip_distance-rho_dip_distance_min(i))/prm%b_sl(i) & - * stt%rho_mob(i,of)*abs(dot_gamma_sl(i)) + * stt%rho_mob(i,me)*abs(dot_gamma_sl(i)) else dot_rho_dip_formation(i) = 0.0_pReal endif @@ -707,39 +707,39 @@ module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Q_cl/(kB*T)) & * (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(kB*T)) - 1.0_pReal) - dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,of) & + dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,me) & / (rho_dip_distance-rho_dip_distance_min(i)) endif endif significantSlipStress enddo slipState - dot%rho_mob(:,of) = abs(dot_gamma_sl)/(prm%b_sl*dst%Lambda_sl(:,of)) & + dot%rho_mob(:,me) = abs(dot_gamma_sl)/(prm%b_sl*dst%Lambda_sl(:,me)) & - dot_rho_dip_formation & - - 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_mob(:,of)*abs(dot_gamma_sl) + - 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_mob(:,me)*abs(dot_gamma_sl) - dot%rho_dip(:,of) = dot_rho_dip_formation & - - 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_dip(:,of)*abs(dot_gamma_sl) & + dot%rho_dip(:,me) = dot_rho_dip_formation & + - 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_dip(:,me)*abs(dot_gamma_sl) & - dot_rho_dip_climb - call kinetics_twin(Mp,T,dot_gamma_sl,instance,of,dot_gamma_twin) - dot%f_tw(:,of) = f_unrotated*dot_gamma_twin/prm%gamma_char + call kinetics_twin(Mp,T,dot_gamma_sl,phase_plasticityInstance(ph),me,dot_gamma_twin) + dot%f_tw(:,me) = f_unrotated*dot_gamma_twin/prm%gamma_char - call kinetics_trans(Mp,T,dot_gamma_sl,instance,of,dot_gamma_tr) - dot%f_tr(:,of) = f_unrotated*dot_gamma_tr + call kinetics_trans(Mp,T,dot_gamma_sl,phase_plasticityInstance(ph),me,dot_gamma_tr) + dot%f_tr(:,me) = f_unrotated*dot_gamma_tr end associate -end subroutine plastic_dislotwin_dotState +end subroutine dislotwin_dotState !-------------------------------------------------------------------------------------------------- !> @brief Calculate derived quantities from state. !-------------------------------------------------------------------------------------------------- -module subroutine plastic_dislotwin_dependentState(T,instance,of) +module subroutine dislotwin_dependentState(T,instance,me) integer, intent(in) :: & instance, & - of + me real(pReal), intent(in) :: & T @@ -763,18 +763,18 @@ module subroutine plastic_dislotwin_dependentState(T,instance,of) stt => state(instance),& dst => dependentState(instance)) - sumf_twin = sum(stt%f_tw(1:prm%sum_N_tw,of)) - sumf_trans = sum(stt%f_tr(1:prm%sum_N_tr,of)) + sumf_twin = sum(stt%f_tw(1:prm%sum_N_tw,me)) + sumf_trans = sum(stt%f_tr(1:prm%sum_N_tr,me)) Gamma = prm%Gamma_sf_0K + prm%dGamma_sf_dT * T !* rescaled volume fraction for topology - f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,of)/prm%t_tw ! this is per system ... + f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,me)/prm%t_tw ! this is per system ... f_over_t_tr = sumf_trans/prm%t_tr ! but this not ! ToDo ...Physically correct, but naming could be adjusted inv_lambda_sl_sl = sqrt(matmul(prm%forestProjection, & - stt%rho_mob(:,of)+stt%rho_dip(:,of)))/prm%i_sl + stt%rho_mob(:,me)+stt%rho_dip(:,me)))/prm%i_sl if (prm%sum_N_tw > 0 .and. prm%sum_N_sl > 0) & inv_lambda_sl_tw = matmul(prm%h_sl_tw,f_over_t_tw)/(1.0_pReal-sumf_twin) @@ -787,41 +787,41 @@ module subroutine plastic_dislotwin_dependentState(T,instance,of) inv_lambda_tr_tr = matmul(prm%h_tr_tr,f_over_t_tr)/(1.0_pReal-sumf_trans) if ((prm%sum_N_tw > 0) .or. (prm%sum_N_tr > 0)) then ! ToDo: better logic needed here - dst%Lambda_sl(:,of) = prm%D & + dst%Lambda_sl(:,me) = prm%D & / (1.0_pReal+prm%D*(inv_lambda_sl_sl + inv_lambda_sl_tw + inv_lambda_sl_tr)) else - dst%Lambda_sl(:,of) = prm%D & + dst%Lambda_sl(:,me) = prm%D & / (1.0_pReal+prm%D*inv_lambda_sl_sl) !!!!!! correct? endif - dst%Lambda_tw(:,of) = prm%i_tw*prm%D/(1.0_pReal+prm%D*inv_lambda_tw_tw) - dst%Lambda_tr(:,of) = prm%i_tr*prm%D/(1.0_pReal+prm%D*inv_lambda_tr_tr) + dst%Lambda_tw(:,me) = prm%i_tw*prm%D/(1.0_pReal+prm%D*inv_lambda_tw_tw) + dst%Lambda_tr(:,me) = prm%i_tr*prm%D/(1.0_pReal+prm%D*inv_lambda_tr_tr) !* threshold stress for dislocation motion - dst%tau_pass(:,of) = prm%mu*prm%b_sl* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,of)+stt%rho_dip(:,of))) + dst%tau_pass(:,me) = prm%mu*prm%b_sl* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,me)+stt%rho_dip(:,me))) !* threshold stress for growing twin/martensite if(prm%sum_N_tw == prm%sum_N_sl) & - dst%tau_hat_tw(:,of) = Gamma/(3.0_pReal*prm%b_tw) & + dst%tau_hat_tw(:,me) = Gamma/(3.0_pReal*prm%b_tw) & + 3.0_pReal*prm%b_tw*prm%mu/(prm%L_tw*prm%b_sl) ! slip Burgers here correct? if(prm%sum_N_tr == prm%sum_N_sl) & - dst%tau_hat_tr(:,of) = Gamma/(3.0_pReal*prm%b_tr) & + dst%tau_hat_tr(:,me) = Gamma/(3.0_pReal*prm%b_tr) & + 3.0_pReal*prm%b_tr*prm%mu/(prm%L_tr*prm%b_sl) & ! slip Burgers here correct? + prm%h*prm%delta_G/ (3.0_pReal*prm%b_tr) - dst%V_tw(:,of) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,of)**2.0_pReal - dst%V_tr(:,of) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,of)**2.0_pReal + dst%V_tw(:,me) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,me)**2.0_pReal + dst%V_tr(:,me) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,me)**2.0_pReal x0 = prm%mu*prm%b_tw**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip and is the same for twin and trans - dst%tau_r_tw(:,of) = prm%mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(pi/3.0_pReal)/x0) + dst%tau_r_tw(:,me) = prm%mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(pi/3.0_pReal)/x0) x0 = prm%mu*prm%b_tr**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip - dst%tau_r_tr(:,of) = prm%mu*prm%b_tr/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(pi/3.0_pReal)/x0) + dst%tau_r_tr(:,me) = prm%mu*prm%b_tr/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(pi/3.0_pReal)/x0) end associate -end subroutine plastic_dislotwin_dependentState +end subroutine dislotwin_dependentState !-------------------------------------------------------------------------------------------------- @@ -882,7 +882,7 @@ end subroutine plastic_dislotwin_results ! 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(Mp,T,instance,of, & +pure subroutine kinetics_slip(Mp,T,instance,me, & dot_gamma_sl,ddot_gamma_dtau_slip,tau_slip) real(pReal), dimension(3,3), intent(in) :: & @@ -891,7 +891,7 @@ pure subroutine kinetics_slip(Mp,T,instance,of, & T !< temperature integer, intent(in) :: & instance, & - of + me real(pReal), dimension(param(instance)%sum_N_sl), intent(out) :: & dot_gamma_sl @@ -920,7 +920,7 @@ pure subroutine kinetics_slip(Mp,T,instance,of, & tau(i) = math_tensordot(Mp,prm%P_sl(1:3,1:3,i)) enddo - tau_eff = abs(tau)-dst%tau_pass(:,of) + tau_eff = abs(tau)-dst%tau_pass(:,me) significantStress: where(tau_eff > tol_math_check) stressRatio = tau_eff/prm%tau_0 @@ -929,7 +929,7 @@ pure subroutine kinetics_slip(Mp,T,instance,of, & v_wait_inverse = prm%v_0**(-1.0_pReal) * exp(BoltzmannRatio*(1.0_pReal-StressRatio_p)** prm%q) v_run_inverse = prm%B/(tau_eff*prm%b_sl) - dot_gamma_sl = sign(stt%rho_mob(:,of)*prm%b_sl/(v_wait_inverse+v_run_inverse),tau) + dot_gamma_sl = sign(stt%rho_mob(:,me)*prm%b_sl/(v_wait_inverse+v_run_inverse),tau) dV_wait_inverse_dTau = -1.0_pReal * v_wait_inverse * prm%p * prm%q * BoltzmannRatio & * (stressRatio**(prm%p-1.0_pReal)) & @@ -938,7 +938,7 @@ pure subroutine kinetics_slip(Mp,T,instance,of, & dV_run_inverse_dTau = -1.0_pReal * v_run_inverse/tau_eff dV_dTau = -1.0_pReal * (dV_wait_inverse_dTau+dV_run_inverse_dTau) & / (v_wait_inverse+v_run_inverse)**2.0_pReal - ddot_gamma_dtau = dV_dTau*stt%rho_mob(:,of)*prm%b_sl + ddot_gamma_dtau = dV_dTau*stt%rho_mob(:,me)*prm%b_sl else where significantStress dot_gamma_sl = 0.0_pReal ddot_gamma_dtau = 0.0_pReal @@ -959,7 +959,7 @@ end subroutine kinetics_slip ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- -pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,of,& +pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,me,& dot_gamma_twin,ddot_gamma_dtau_twin) real(pReal), dimension(3,3), intent(in) :: & @@ -968,7 +968,7 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,of,& T !< temperature integer, intent(in) :: & instance, & - of + me real(pReal), dimension(param(instance)%sum_N_sl), intent(in) :: & dot_gamma_sl @@ -992,11 +992,11 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,of,& isFCC: if (prm%fccTwinTransNucleation) then s1=prm%fcc_twinNucleationSlipPair(1,i) s2=prm%fcc_twinNucleationSlipPair(2,i) - if (tau(i) < dst%tau_r_tw(i,of)) then ! ToDo: correct? - Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,of)+stt%rho_dip(s2,of))+& - abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,of)+stt%rho_dip(s1,of)))/& ! ToDo: MD: it would be more consistent to use shearrates from state + if (tau(i) < dst%tau_r_tw(i,me)) then ! ToDo: correct? + Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,me)+stt%rho_dip(s2,me))+& + abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,me)+stt%rho_dip(s1,me)))/& ! ToDo: MD: it would be more consistent to use shearrates from state (prm%L_tw*prm%b_sl(i))*& - (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tw(i,of)-tau(i)))) ! P_ncs + (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tw(i,me)-tau(i)))) ! P_ncs else Ndot0=0.0_pReal end if @@ -1006,8 +1006,8 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,of,& enddo significantStress: where(tau > tol_math_check) - StressRatio_r = (dst%tau_hat_tw(:,of)/tau)**prm%r - dot_gamma_twin = prm%gamma_char * dst%V_tw(:,of) * Ndot0*exp(-StressRatio_r) + StressRatio_r = (dst%tau_hat_tw(:,me)/tau)**prm%r + dot_gamma_twin = prm%gamma_char * dst%V_tw(:,me) * Ndot0*exp(-StressRatio_r) ddot_gamma_dtau = (dot_gamma_twin*prm%r/tau)*StressRatio_r else where significantStress dot_gamma_twin = 0.0_pReal @@ -1028,7 +1028,7 @@ end subroutine kinetics_twin ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- -pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,of,& +pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,me,& dot_gamma_tr,ddot_gamma_dtau_trans) real(pReal), dimension(3,3), intent(in) :: & @@ -1037,7 +1037,7 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,of,& T !< temperature integer, intent(in) :: & instance, & - of + me real(pReal), dimension(param(instance)%sum_N_sl), intent(in) :: & dot_gamma_sl @@ -1060,11 +1060,11 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,of,& isFCC: if (prm%fccTwinTransNucleation) then s1=prm%fcc_twinNucleationSlipPair(1,i) s2=prm%fcc_twinNucleationSlipPair(2,i) - if (tau(i) < dst%tau_r_tr(i,of)) then ! ToDo: correct? - Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,of)+stt%rho_dip(s2,of))+& - abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,of)+stt%rho_dip(s1,of)))/& ! ToDo: MD: it would be more consistent to use shearrates from state + if (tau(i) < dst%tau_r_tr(i,me)) then ! ToDo: correct? + Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,me)+stt%rho_dip(s2,me))+& + abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,me)+stt%rho_dip(s1,me)))/& ! ToDo: MD: it would be more consistent to use shearrates from state (prm%L_tr*prm%b_sl(i))*& - (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tr(i,of)-tau(i)))) ! P_ncs + (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tr(i,me)-tau(i)))) ! P_ncs else Ndot0=0.0_pReal end if @@ -1074,8 +1074,8 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,of,& enddo significantStress: where(tau > tol_math_check) - StressRatio_s = (dst%tau_hat_tr(:,of)/tau)**prm%s - dot_gamma_tr = dst%V_tr(:,of) * Ndot0*exp(-StressRatio_s) + StressRatio_s = (dst%tau_hat_tr(:,me)/tau)**prm%s + dot_gamma_tr = dst%V_tr(:,me) * Ndot0*exp(-StressRatio_s) ddot_gamma_dtau = (dot_gamma_tr*prm%s/tau)*StressRatio_s else where significantStress dot_gamma_tr = 0.0_pReal @@ -1088,4 +1088,4 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,of,& end subroutine kinetics_trans -end submodule plastic_dislotwin +end submodule dislotwin diff --git a/src/constitutive_plastic_isotropic.f90 b/src/phase_mechanics_plastic_isotropic.f90 similarity index 86% rename from src/constitutive_plastic_isotropic.f90 rename to src/phase_mechanics_plastic_isotropic.f90 index b7c5f67c1..6789e74b4 100644 --- a/src/constitutive_plastic_isotropic.f90 +++ b/src/phase_mechanics_plastic_isotropic.f90 @@ -7,7 +7,7 @@ !! resolving the stress on the slip systems. Will give the response of phenopowerlaw for an !! untextured polycrystal !-------------------------------------------------------------------------------------------------- -submodule(constitutive:constitutive_mech) plastic_isotropic +submodule(phase:plastic) isotropic type :: tParameters real(pReal) :: & @@ -68,7 +68,7 @@ module function plastic_isotropic_init() result(myPlasticity) mech, & pl - print'(/,a)', ' <<<+- plastic_isotropic init -+>>>' + print'(/,a)', ' <<<+- phase:mechanics:plastic:isotropic init -+>>>' myPlasticity = plastic_active('isotropic') Ninstances = count(myPlasticity) @@ -131,7 +131,7 @@ module function plastic_isotropic_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! allocate state arrays - Nconstituents = count(material_phaseAt == p) * discretization_nIPs + Nconstituents = count(material_phaseAt2 == p) sizeDotState = size(['xi ','gamma']) sizeState = sizeDotState @@ -168,7 +168,7 @@ end function plastic_isotropic_init !-------------------------------------------------------------------------------------------------- !> @brief Calculate plastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- -module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) +module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient @@ -178,8 +178,8 @@ module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & - instance, & - of + ph, & + me real(pReal), dimension(3,3) :: & Mp_dev !< deviatoric part of the Mandel stress @@ -190,24 +190,16 @@ module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) integer :: & k, l, m, n - associate(prm => param(instance), stt => state(instance)) + associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph))) Mp_dev = math_deviatoric33(Mp) squarenorm_Mp_dev = math_tensordot(Mp_dev,Mp_dev) norm_Mp_dev = sqrt(squarenorm_Mp_dev) if (norm_Mp_dev > 0.0_pReal) then - dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%M*stt%xi(of))) **prm%n + dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%M*stt%xi(me))) **prm%n Lp = dot_gamma/prm%M * Mp_dev/norm_Mp_dev -#ifdef DEBUG - if (debugConstitutive%extensive .and. (of == prm%of_debug .or. .not. debugConstitutive%selective)) then - print'(/,a,/,3(12x,3(f12.4,1x)/))', '<< CONST isotropic >> Tstar (dev) / MPa', & - transpose(Mp_dev)*1.0e-6_pReal - print'(/,a,/,f12.5)', '<< CONST isotropic >> norm Tstar / MPa', norm_Mp_dev*1.0e-6_pReal - print'(/,a,/,f12.5)', '<< CONST isotropic >> gdot', dot_gamma - end if -#endif forall (k=1:3,l=1:3,m=1:3,n=1:3) & dLp_dMp(k,l,m,n) = (prm%n-1.0_pReal) * Mp_dev(k,l)*Mp_dev(m,n) / squarenorm_Mp_dev forall (k=1:3,l=1:3) & @@ -222,13 +214,13 @@ module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) end associate -end subroutine plastic_isotropic_LpAndItsTangent +end subroutine isotropic_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief Calculate inelastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- -module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of) +module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,me) real(pReal), dimension(3,3), intent(out) :: & Li !< inleastic velocity gradient @@ -239,7 +231,7 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of) Mi !< Mandel stress integer, intent(in) :: & instance, & - of + me real(pReal) :: & tr !< trace of spherical part of Mandel stress (= 3 x pressure) @@ -252,19 +244,10 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of) if (prm%dilatation .and. abs(tr) > 0.0_pReal) then ! no stress or J2 plasticity --> Li and its derivative are zero Li = math_I3 & - * prm%dot_gamma_0/prm%M * (3.0_pReal*prm%M*stt%xi(of))**(-prm%n) & + * prm%dot_gamma_0/prm%M * (3.0_pReal*prm%M*stt%xi(me))**(-prm%n) & * tr * abs(tr)**(prm%n-1.0_pReal) -#ifdef DEBUG - if (debugConstitutive%extensive .and. (of == prm%of_debug .or. .not. debugConstitutive%selective)) then - print'(/,a,/,f12.5)', '<< CONST isotropic >> pressure / MPa', tr/3.0_pReal*1.0e-6_pReal - print'(/,a,/,f12.5)', '<< CONST isotropic >> gdot', prm%dot_gamma_0 * (3.0_pReal*prm%M*stt%xi(of))**(-prm%n) & - * tr * abs(tr)**(prm%n-1.0_pReal) - end if -#endif - - forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLi_dMi(k,l,m,n) = prm%n / tr * Li(k,l) * math_I3(m,n) + forall (k=1:3,l=1:3,m=1:3,n=1:3) dLi_dMi(k,l,m,n) = prm%n / tr * Li(k,l) * math_I3(m,n) else Li = 0.0_pReal @@ -279,20 +262,21 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of) !-------------------------------------------------------------------------------------------------- !> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- -module subroutine plastic_isotropic_dotState(Mp,instance,of) +module subroutine isotropic_dotState(Mp,ph,me) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & - instance, & - of + ph, & + me real(pReal) :: & dot_gamma, & !< strainrate xi_inf_star, & !< saturation xi norm_Mp !< norm of the (deviatoric) Mandel stress - associate(prm => param(instance), stt => state(instance), dot => dotState(instance)) + associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph)), & + dot => dotState(phase_plasticityInstance(ph))) if (prm%dilatation) then norm_Mp = sqrt(math_tensordot(Mp,Mp)) @@ -300,7 +284,7 @@ module subroutine plastic_isotropic_dotState(Mp,instance,of) norm_Mp = sqrt(math_tensordot(math_deviatoric33(Mp),math_deviatoric33(Mp))) endif - dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp /(prm%M*stt%xi(of))) **prm%n + dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp /(prm%M*stt%xi(me))) **prm%n if (dot_gamma > 1e-12_pReal) then if (dEq0(prm%c_1)) then @@ -310,19 +294,19 @@ module subroutine plastic_isotropic_dotState(Mp,instance,of) + asinh( (dot_gamma / prm%c_1)**(1.0_pReal / prm%c_2))**(1.0_pReal / prm%c_3) & / prm%c_4 * (dot_gamma / prm%dot_gamma_0)**(1.0_pReal / prm%n) endif - dot%xi(of) = dot_gamma & + dot%xi(me) = dot_gamma & * ( prm%h_0 + prm%h_ln * log(dot_gamma) ) & - * abs( 1.0_pReal - stt%xi(of)/xi_inf_star )**prm%a & - * sign(1.0_pReal, 1.0_pReal - stt%xi(of)/xi_inf_star) + * abs( 1.0_pReal - stt%xi(me)/xi_inf_star )**prm%a & + * sign(1.0_pReal, 1.0_pReal - stt%xi(me)/xi_inf_star) else - dot%xi(of) = 0.0_pReal + dot%xi(me) = 0.0_pReal endif - dot%gamma(of) = dot_gamma ! ToDo: not really used + dot%gamma(me) = dot_gamma ! ToDo: not really used end associate -end subroutine plastic_isotropic_dotState +end subroutine isotropic_dotState !-------------------------------------------------------------------------------------------------- @@ -348,4 +332,4 @@ module subroutine plastic_isotropic_results(instance,group) end subroutine plastic_isotropic_results -end submodule plastic_isotropic +end submodule isotropic diff --git a/src/constitutive_plastic_kinehardening.f90 b/src/phase_mechanics_plastic_kinehardening.f90 similarity index 89% rename from src/constitutive_plastic_kinehardening.f90 rename to src/phase_mechanics_plastic_kinehardening.f90 index 8454b28f8..919302bc1 100644 --- a/src/constitutive_plastic_kinehardening.f90 +++ b/src/phase_mechanics_plastic_kinehardening.f90 @@ -5,7 +5,7 @@ !> @brief Phenomenological crystal plasticity using a power law formulation for the shear rates !! and a Voce-type kinematic hardening rule !-------------------------------------------------------------------------------------------------- -submodule(constitutive:constitutive_mech) plastic_kinehardening +submodule(phase:plastic) kinehardening type :: tParameters real(pReal) :: & @@ -80,7 +80,7 @@ module function plastic_kinehardening_init() result(myPlasticity) mech, & pl - print'(/,a)', ' <<<+- plastic_kinehardening init -+>>>' + print'(/,a)', ' <<<+- phase:mechanics:plastic:kinehardening init -+>>>' myPlasticity = plastic_active('kinehardening') Ninstances = count(myPlasticity) @@ -175,7 +175,7 @@ module function plastic_kinehardening_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! allocate state arrays - Nconstituents = count(material_phaseAt == p) * discretization_nIPs + Nconstituents = count(material_phaseAt2 == p) sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%sum_N_sl!ToDo: adjust names, ask Philip sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl !ToDo: adjust names sizeState = sizeDotState + sizeDeltaState @@ -240,7 +240,7 @@ end function plastic_kinehardening_init !-------------------------------------------------------------------------------------------------- !> @brief Calculate plastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- -pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) +pure module subroutine kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient @@ -250,21 +250,21 @@ pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,insta real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & - instance, & - of + ph, & + me integer :: & i,k,l,m,n - real(pReal), dimension(param(instance)%sum_N_sl) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & gdot_pos,gdot_neg, & dgdot_dtau_pos,dgdot_dtau_neg Lp = 0.0_pReal dLp_dMp = 0.0_pReal - associate(prm => param(instance)) + associate(prm => param(phase_plasticityInstance(ph))) - call kinetics(Mp,instance,of,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) + call kinetics(Mp,phase_plasticityInstance(ph),me,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) do i = 1, prm%sum_N_sl Lp = Lp + (gdot_pos(i)+gdot_neg(i))*prm%P(1:3,1:3,i) @@ -276,44 +276,45 @@ pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,insta end associate -end subroutine plastic_kinehardening_LpAndItsTangent +end subroutine kinehardening_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- -module subroutine plastic_kinehardening_dotState(Mp,instance,of) +module subroutine plastic_kinehardening_dotState(Mp,ph,me) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & - instance, & - of + ph, & + me real(pReal) :: & sumGamma - real(pReal), dimension(param(instance)%sum_N_sl) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & gdot_pos,gdot_neg - associate(prm => param(instance), stt => state(instance), dot => dotState(instance)) + associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph)),& + dot => dotState(phase_plasticityInstance(ph))) - call kinetics(Mp,instance,of,gdot_pos,gdot_neg) - dot%accshear(:,of) = abs(gdot_pos+gdot_neg) - sumGamma = sum(stt%accshear(:,of)) + call kinetics(Mp,phase_plasticityInstance(ph),me,gdot_pos,gdot_neg) + dot%accshear(:,me) = abs(gdot_pos+gdot_neg) + sumGamma = sum(stt%accshear(:,me)) - dot%crss(:,of) = matmul(prm%interaction_SlipSlip,dot%accshear(:,of)) & + dot%crss(:,me) = matmul(prm%interaction_SlipSlip,dot%accshear(:,me)) & * ( prm%h_inf_f & + (prm%h_0_f - prm%h_inf_f + prm%h_0_f*prm%h_inf_f*sumGamma/prm%xi_inf_f) & * exp(-sumGamma*prm%h_0_f/prm%xi_inf_f) & ) - dot%crss_back(:,of) = stt%sense(:,of)*dot%accshear(:,of) * & + dot%crss_back(:,me) = stt%sense(:,me)*dot%accshear(:,me) * & ( prm%h_inf_b + & (prm%h_0_b - prm%h_inf_b & - + prm%h_0_b*prm%h_inf_b/(prm%xi_inf_b+stt%chi0(:,of))*(stt%accshear(:,of)-stt%gamma0(:,of))& - ) *exp(-(stt%accshear(:,of)-stt%gamma0(:,of)) *prm%h_0_b/(prm%xi_inf_b+stt%chi0(:,of))) & + + prm%h_0_b*prm%h_inf_b/(prm%xi_inf_b+stt%chi0(:,me))*(stt%accshear(:,me)-stt%gamma0(:,me))& + ) *exp(-(stt%accshear(:,me)-stt%gamma0(:,me)) *prm%h_0_b/(prm%xi_inf_b+stt%chi0(:,me))) & ) end associate @@ -324,13 +325,13 @@ end subroutine plastic_kinehardening_dotState !-------------------------------------------------------------------------------------------------- !> @brief Calculate (instantaneous) incremental change of microstructure. !-------------------------------------------------------------------------------------------------- -module subroutine plastic_kinehardening_deltaState(Mp,instance,of) +module subroutine plastic_kinehardening_deltaState(Mp,instance,me) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & instance, & - of + me real(pReal), dimension(param(instance)%sum_N_sl) :: & gdot_pos,gdot_neg, & @@ -338,29 +339,29 @@ module subroutine plastic_kinehardening_deltaState(Mp,instance,of) associate(prm => param(instance), stt => state(instance), dlt => deltaState(instance)) - call kinetics(Mp,instance,of,gdot_pos,gdot_neg) - sense = merge(state(instance)%sense(:,of), & ! keep existing... + call kinetics(Mp,instance,me,gdot_pos,gdot_neg) + sense = merge(state(instance)%sense(:,me), & ! keep existing... sign(1.0_pReal,gdot_pos+gdot_neg), & ! ...or have a defined dEq0(gdot_pos+gdot_neg,1e-10_pReal)) ! current sense of shear direction #ifdef DEBUG if (debugConstitutive%extensive & - .and. (of == prm%of_debug .or. .not. debugConstitutive%selective)) then + .and. (me == prm%of_debug .or. .not. debugConstitutive%selective)) then print*, '======= kinehardening delta state =======' - print*, sense,state(instance)%sense(:,of) + print*, sense,state(instance)%sense(:,me) endif #endif !-------------------------------------------------------------------------------------------------- -! switch in sense of shear? - where(dNeq(sense,stt%sense(:,of),0.1_pReal)) - dlt%sense (:,of) = sense - stt%sense(:,of) ! switch sense - dlt%chi0 (:,of) = abs(stt%crss_back(:,of)) - stt%chi0(:,of) ! remember current backstress magnitude - dlt%gamma0(:,of) = stt%accshear(:,of) - stt%gamma0(:,of) ! remember current accumulated shear +! switch in sense me shear? + where(dNeq(sense,stt%sense(:,me),0.1_pReal)) + dlt%sense (:,me) = sense - stt%sense(:,me) ! switch sense + dlt%chi0 (:,me) = abs(stt%crss_back(:,me)) - stt%chi0(:,me) ! remember current backstress magnitude + dlt%gamma0(:,me) = stt%accshear(:,me) - stt%gamma0(:,me) ! remember current accumulated shear else where - dlt%sense (:,of) = 0.0_pReal - dlt%chi0 (:,of) = 0.0_pReal - dlt%gamma0(:,of) = 0.0_pReal + dlt%sense (:,me) = 0.0_pReal + dlt%chi0 (:,me) = 0.0_pReal + dlt%gamma0(:,me) = 0.0_pReal end where end associate @@ -413,14 +414,14 @@ end subroutine plastic_kinehardening_results ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- -pure subroutine kinetics(Mp,instance,of, & +pure subroutine kinetics(Mp,instance,me, & gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & instance, & - of + me real(pReal), intent(out), dimension(param(instance)%sum_N_sl) :: & gdot_pos, & @@ -437,21 +438,21 @@ pure subroutine kinetics(Mp,instance,of, & associate(prm => param(instance), stt => state(instance)) do i = 1, prm%sum_N_sl - tau_pos(i) = math_tensordot(Mp,prm%nonSchmid_pos(1:3,1:3,i)) - stt%crss_back(i,of) - tau_neg(i) = merge(math_tensordot(Mp,prm%nonSchmid_neg(1:3,1:3,i)) - stt%crss_back(i,of), & + tau_pos(i) = math_tensordot(Mp,prm%nonSchmid_pos(1:3,1:3,i)) - stt%crss_back(i,me) + tau_neg(i) = merge(math_tensordot(Mp,prm%nonSchmid_neg(1:3,1:3,i)) - stt%crss_back(i,me), & 0.0_pReal, prm%nonSchmidActive) enddo where(dNeq0(tau_pos)) gdot_pos = prm%dot_gamma_0 * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active - * sign(abs(tau_pos/stt%crss(:,of))**prm%n, tau_pos) + * sign(abs(tau_pos/stt%crss(:,me))**prm%n, tau_pos) else where gdot_pos = 0.0_pReal end where where(dNeq0(tau_neg)) gdot_neg = prm%dot_gamma_0 * 0.5_pReal & ! only used if non-Schmid active, always 1/2 - * sign(abs(tau_neg/stt%crss(:,of))**prm%n, tau_neg) + * sign(abs(tau_neg/stt%crss(:,me))**prm%n, tau_neg) else where gdot_neg = 0.0_pReal end where @@ -474,4 +475,4 @@ pure subroutine kinetics(Mp,instance,of, & end subroutine kinetics -end submodule plastic_kinehardening +end submodule kinehardening diff --git a/src/constitutive_plastic_none.f90 b/src/phase_mechanics_plastic_none.f90 similarity index 89% rename from src/constitutive_plastic_none.f90 rename to src/phase_mechanics_plastic_none.f90 index 27a01fb93..b8d1678eb 100644 --- a/src/constitutive_plastic_none.f90 +++ b/src/phase_mechanics_plastic_none.f90 @@ -4,7 +4,7 @@ !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief Dummy plasticity for purely elastic material !-------------------------------------------------------------------------------------------------- -submodule(constitutive:constitutive_mech) plastic_none +submodule(phase:plastic) none contains @@ -25,7 +25,7 @@ module function plastic_none_init() result(myPlasticity) mech, & pl - print'(/,a)', ' <<<+- plastic_none init -+>>>' + print'(/,a)', ' <<<+- phase:mechanics:plastic:none init -+>>>' phases => config_material%get('phase') allocate(myPlasticity(phases%length), source = .false.) @@ -43,11 +43,11 @@ module function plastic_none_init() result(myPlasticity) do p = 1, phases%length phase => phases%get(p) if(.not. myPlasticity(p)) cycle - Nconstituents = count(material_phaseAt == p) * discretization_nIPs + Nconstituents = count(material_phaseAt2 == p) call constitutive_allocateState(plasticState(p),Nconstituents,0,0,0) enddo end function plastic_none_init -end submodule plastic_none +end submodule none diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/phase_mechanics_plastic_nonlocal.f90 similarity index 96% rename from src/constitutive_plastic_nonlocal.f90 rename to src/phase_mechanics_plastic_nonlocal.f90 index 2244eb7ad..693ffcf93 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/phase_mechanics_plastic_nonlocal.f90 @@ -4,7 +4,7 @@ !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @brief material subroutine for plasticity including dislocation flux !-------------------------------------------------------------------------------------------------- -submodule(constitutive:constitutive_mech) plastic_nonlocal +submodule(phase:plastic) nonlocal use geometry_plastic_nonlocal, only: & nIPneighbors => geometry_plastic_nonlocal_nIPneighbors, & IPneighborhood => geometry_plastic_nonlocal_IPneighborhood, & @@ -187,7 +187,7 @@ module function plastic_nonlocal_init() result(myPlasticity) mech, & pl - print'(/,a)', ' <<<+- plastic_nonlocal init -+>>>' + print'(/,a)', ' <<<+- phase:mechanics:plastic:nonlocal init -+>>>' myPlasticity = plastic_active('nonlocal') Ninstances = count(myPlasticity) @@ -393,7 +393,7 @@ module function plastic_nonlocal_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! allocate state arrays - Nconstituents = count(material_phaseAt==p) * discretization_nIPs + Nconstituents = count(material_phaseAt2 == p) sizeDotState = size([ 'rhoSglEdgePosMobile ','rhoSglEdgeNegMobile ', & 'rhoSglScrewPosMobile ','rhoSglScrewNegMobile ', & 'rhoSglEdgePosImmobile ','rhoSglEdgeNegImmobile ', & @@ -552,17 +552,16 @@ end function plastic_nonlocal_init !-------------------------------------------------------------------------------------------------- !> @brief calculates quantities characterizing the microstructure !-------------------------------------------------------------------------------------------------- -module subroutine plastic_nonlocal_dependentState(instance, of, ip, el) +module subroutine nonlocal_dependentState(instance, me, ip, el) integer, intent(in) :: & instance, & - of, & + me, & ip, & el integer :: & ph, & - me, & no, & !< neighbor offset neighbor_el, & ! element number of neighboring material point neighbor_ip, & ! integration point of neighboring material point @@ -612,9 +611,9 @@ module subroutine plastic_nonlocal_dependentState(instance, of, ip, el) associate(prm => param(instance),dst => microstructure(instance), stt => state(instance)) - rho = getRho(instance,of,ip,el) + rho = getRho(instance,me,ip,el) - stt%rho_forest(:,of) = matmul(prm%forestProjection_Edge, sum(abs(rho(:,edg)),2)) & + stt%rho_forest(:,me) = matmul(prm%forestProjection_Edge, sum(abs(rho(:,edg)),2)) & + matmul(prm%forestProjection_Screw,sum(abs(rho(:,scr)),2)) @@ -624,13 +623,13 @@ module subroutine plastic_nonlocal_dependentState(instance, of, ip, el) myInteractionMatrix = prm%h_sl_sl & * spread(( 1.0_pReal - prm%f_F & + prm%f_F & - * log(0.35_pReal * prm%b_sl * sqrt(max(stt%rho_forest(:,of),prm%rho_significant))) & + * log(0.35_pReal * prm%b_sl * sqrt(max(stt%rho_forest(:,me),prm%rho_significant))) & / log(0.35_pReal * prm%b_sl * 1e6_pReal))** 2.0_pReal,2,prm%sum_N_sl) else myInteractionMatrix = prm%h_sl_sl endif - dst%tau_pass(:,of) = prm%mu * prm%b_sl & + dst%tau_pass(:,me) = prm%mu * prm%b_sl & * sqrt(matmul(myInteractionMatrix,sum(abs(rho),2))) !*** calculate the dislocation stress of the neighboring excess dislocation densities @@ -640,10 +639,9 @@ module subroutine plastic_nonlocal_dependentState(instance, of, ip, el) ! ToDo: MD: this is most likely only correct for F_i = I !################################################################################################# - rho0 = getRho0(instance,of,ip,el) + rho0 = getRho0(instance,me,ip,el) if (.not. phase_localPlasticity(material_phaseAt(1,el)) .and. prm%shortRangeStressCorrection) then ph = material_phaseAt(1,el) - me = material_phaseMemberAt(1,ip,el) invFp = math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,me)) invFe = math_inv33(constitutive_mech_Fe(ph)%data(1:3,1:3,me)) @@ -734,7 +732,7 @@ module subroutine plastic_nonlocal_dependentState(instance, of, ip, el) where(rhoTotal > 0.0_pReal) rhoExcessGradient_over_rho = rhoExcessGradient / rhoTotal ! ... gives the local stress correction when multiplied with a factor - dst%tau_back(s,of) = - prm%mu * prm%b_sl(s) / (2.0_pReal * PI) & + dst%tau_back(s,me) = - prm%mu * prm%b_sl(s) / (2.0_pReal * PI) & * ( rhoExcessGradient_over_rho(1) / (1.0_pReal - prm%nu) & + rhoExcessGradient_over_rho(2)) enddo @@ -745,29 +743,29 @@ module subroutine plastic_nonlocal_dependentState(instance, of, ip, el) .and. ((debugConstitutive%element == el .and. debugConstitutive%ip == ip)& .or. .not. debugConstitutive%selective)) then print'(/,a,i8,1x,i2,1x,i1,/)', '<< CONST >> nonlocal_microstructure at el ip ',el,ip - print'(a,/,12x,12(e10.3,1x))', '<< CONST >> rhoForest', stt%rho_forest(:,of) - print'(a,/,12x,12(f10.5,1x))', '<< CONST >> tauThreshold / MPa', dst%tau_pass(:,of)*1e-6 - print'(a,/,12x,12(f10.5,1x),/)', '<< CONST >> tauBack / MPa', dst%tau_back(:,of)*1e-6 + print'(a,/,12x,12(e10.3,1x))', '<< CONST >> rhoForest', stt%rho_forest(:,me) + print'(a,/,12x,12(f10.5,1x))', '<< CONST >> tauThreshold / MPa', dst%tau_pass(:,me)*1e-6 + print'(a,/,12x,12(f10.5,1x),/)', '<< CONST >> tauBack / MPa', dst%tau_back(:,me)*1e-6 endif #endif end associate -end subroutine plastic_nonlocal_dependentState +end subroutine nonlocal_dependentState !-------------------------------------------------------------------------------------------------- !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- -module subroutine plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp, & - Mp,Temperature,instance,of,ip,el) +module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & + Mp,Temperature,ph,me,ip,el) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp integer, intent(in) :: & - instance, & - of, & + ph, & + me, & ip, & !< current integration point el !< current element number real(pReal), intent(in) :: & @@ -784,24 +782,25 @@ module subroutine plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp, & l, & t, & !< dislocation type s !< index of my current slip system - real(pReal), dimension(param(instance)%sum_N_sl,8) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,8) :: & rhoSgl !< single dislocation densities (including blocked) - real(pReal), dimension(param(instance)%sum_N_sl,10) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,10) :: & rho - real(pReal), dimension(param(instance)%sum_N_sl,4) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,4) :: & v, & !< velocity tauNS, & !< resolved shear stress including non Schmid and backstress terms dv_dtau, & !< velocity derivative with respect to the shear stress dv_dtauNS !< velocity derivative with respect to the shear stress - real(pReal), dimension(param(instance)%sum_N_sl) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & tau, & !< resolved shear stress including backstress terms gdotTotal !< shear rate - associate(prm => param(instance),dst=>microstructure(instance),stt=>state(instance)) + associate(prm => param(phase_plasticityInstance(ph)),dst=>microstructure(phase_plasticityInstance(ph)),& + stt=>state(phase_plasticityInstance(ph))) ns = prm%sum_N_sl !*** shortcut to state variables - rho = getRho(instance,of,ip,el) + rho = getRho(phase_plasticityInstance(ph),me,ip,el) rhoSgl = rho(:,sgl) do s = 1,ns @@ -816,12 +815,12 @@ module subroutine plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp, & tauNS(s,4) = math_tensordot(Mp, -prm%nonSchmid_pos(1:3,1:3,s)) endif enddo - tauNS = tauNS + spread(dst%tau_back(:,of),2,4) - tau = tau + dst%tau_back(:,of) + tauNS = tauNS + spread(dst%tau_back(:,me),2,4) + tau = tau + dst%tau_back(:,me) ! edges call kinetics(v(:,1), dv_dtau(:,1), dv_dtauNS(:,1), & - tau, tauNS(:,1), dst%tau_pass(:,of),1,Temperature, instance) + tau, tauNS(:,1), dst%tau_pass(:,me),1,Temperature, phase_plasticityInstance(ph)) v(:,2) = v(:,1) dv_dtau(:,2) = dv_dtau(:,1) dv_dtauNS(:,2) = dv_dtauNS(:,1) @@ -834,11 +833,11 @@ module subroutine plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp, & else do t = 3,4 call kinetics(v(:,t), dv_dtau(:,t), dv_dtauNS(:,t), & - tau, tauNS(:,t), dst%tau_pass(:,of),2,Temperature, instance) + tau, tauNS(:,t), dst%tau_pass(:,me),2,Temperature, phase_plasticityInstance(ph)) enddo endif - stt%v(:,of) = pack(v,.true.) + stt%v(:,me) = pack(v,.true.) !*** Bauschinger effect forall (s = 1:ns, t = 5:8, rhoSgl(s,t) * v(s,t-4) < 0.0_pReal) & @@ -861,19 +860,19 @@ module subroutine plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp, & end associate -end subroutine plastic_nonlocal_LpAndItsTangent +end subroutine nonlocal_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief (instantaneous) incremental change of microstructure !-------------------------------------------------------------------------------------------------- -module subroutine plastic_nonlocal_deltaState(Mp,instance,of,ip,el) +module subroutine plastic_nonlocal_deltaState(Mp,instance,me,ip,el) real(pReal), dimension(3,3), intent(in) :: & Mp !< MandelStress integer, intent(in) :: & instance, & ! current instance of this plasticity - of, & !< offset + me, & !< offset ip, & el @@ -904,10 +903,10 @@ module subroutine plastic_nonlocal_deltaState(Mp,instance,of,ip,el) ns = prm%sum_N_sl !*** shortcut to state variables - forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,instance),of) - forall (s = 1:ns, c = 1:2) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,instance),of) + forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,instance),me) + forall (s = 1:ns, c = 1:2) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,instance),me) - rho = getRho(instance,of,ip,el) + rho = getRho(instance,me,ip,el) rhoDip = rho(:,dip) !**************************************************************************** @@ -928,7 +927,7 @@ module subroutine plastic_nonlocal_deltaState(Mp,instance,of,ip,el) !*** calculate limits for stable dipole height do s = 1,prm%sum_N_sl - tau(s) = math_tensordot(Mp, prm%Schmid(1:3,1:3,s)) +dst%tau_back(s,of) + tau(s) = math_tensordot(Mp, prm%Schmid(1:3,1:3,s)) +dst%tau_back(s,me) if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal enddo @@ -952,10 +951,10 @@ module subroutine plastic_nonlocal_deltaState(Mp,instance,of,ip,el) / (dUpperOld(s,c) - prm%minDipoleHeight(s,c)) forall (t=1:4) deltaRhoDipole2SingleStress(:,t) = -0.5_pReal * deltaRhoDipole2SingleStress(:,(t-1)/2+9) - forall (s = 1:ns, c = 1:2) plasticState(ph)%state(iD(s,c,instance),of) = dUpper(s,c) + forall (s = 1:ns, c = 1:2) plasticState(ph)%state(iD(s,c,instance),me) = dUpper(s,c) - plasticState(ph)%deltaState(:,of) = 0.0_pReal - del%rho(:,of) = reshape(deltaRhoRemobilization + deltaRhoDipole2SingleStress, [10*ns]) + plasticState(ph)%deltaState(:,me) = 0.0_pReal + del%rho(:,me) = reshape(deltaRhoRemobilization + deltaRhoDipole2SingleStress, [10*ns]) #ifdef DEBUG if (debugConstitutive%extensive & @@ -974,8 +973,8 @@ end subroutine plastic_nonlocal_deltaState !--------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !--------------------------------------------------------------------------------------------------- -module subroutine plastic_nonlocal_dotState(Mp, Temperature,timestep, & - instance,of,ip,el) +module subroutine nonlocal_dotState(Mp, Temperature,timestep, & + ph,me,ip,el) real(pReal), dimension(3,3), intent(in) :: & Mp !< MandelStress @@ -983,18 +982,17 @@ module subroutine plastic_nonlocal_dotState(Mp, Temperature,timestep, & Temperature, & !< temperature timestep !< substepped crystallite time increment integer, intent(in) :: & - instance, & - of, & + ph, & + me, & ip, & !< current integration point el !< current element number integer :: & - ph, & ns, & !< short notation for the total number of active slip systems c, & !< character of dislocation t, & !< type of dislocation s !< index of my current slip system - real(pReal), dimension(param(instance)%sum_N_sl,10) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,10) :: & rho, & rho0, & !< dislocation density at beginning of time step rhoDot, & !< density evolution @@ -1002,45 +1000,44 @@ module subroutine plastic_nonlocal_dotState(Mp, Temperature,timestep, & rhoDotSingle2DipoleGlide, & !< density evolution by dipole formation (by glide) rhoDotAthermalAnnihilation, & !< density evolution by athermal annihilation rhoDotThermalAnnihilation !< density evolution by thermal annihilation - real(pReal), dimension(param(instance)%sum_N_sl,8) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,8) :: & rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles) my_rhoSgl0 !< single dislocation densities of central ip (positive/negative screw and edge without dipoles) - real(pReal), dimension(param(instance)%sum_N_sl,4) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,4) :: & v, & !< current dislocation glide velocity v0, & gdot !< shear rates - real(pReal), dimension(param(instance)%sum_N_sl) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & tau, & !< current resolved shear stress vClimb !< climb velocity of edge dipoles - real(pReal), dimension(param(instance)%sum_N_sl,2) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,2) :: & rhoDip, & !< current dipole dislocation densities (screw and edge dipoles) dLower, & !< minimum stable dipole distance for edges and screws dUpper !< current maximum stable dipole distance for edges and screws real(pReal) :: & selfDiffusion !< self diffusion - ph = material_phaseAt(1,el) if (timestep <= 0.0_pReal) then plasticState(ph)%dotState = 0.0_pReal return endif - associate(prm => param(instance), & - dst => microstructure(instance), & - dot => dotState(instance), & - stt => state(instance)) + associate(prm => param(phase_plasticityInstance(ph)), & + dst => microstructure(phase_plasticityInstance(ph)), & + dot => dotState(phase_plasticityInstance(ph)), & + stt => state(phase_plasticityInstance(ph))) ns = prm%sum_N_sl tau = 0.0_pReal gdot = 0.0_pReal - rho = getRho(instance,of,ip,el) + rho = getRho(phase_plasticityInstance(ph),me,ip,el) rhoSgl = rho(:,sgl) rhoDip = rho(:,dip) - rho0 = getRho0(instance,of,ip,el) + rho0 = getRho0(phase_plasticityInstance(ph),me,ip,el) my_rhoSgl0 = rho0(:,sgl) - forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,instance),of) + forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,phase_plasticityInstance(ph)),me) gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4) #ifdef DEBUG @@ -1055,7 +1052,7 @@ module subroutine plastic_nonlocal_dotState(Mp, Temperature,timestep, & !**************************************************************************** !*** limits for stable dipole height do s = 1,ns - tau(s) = math_tensordot(Mp, prm%Schmid(1:3,1:3,s)) + dst%tau_back(s,of) + tau(s) = math_tensordot(Mp, prm%Schmid(1:3,1:3,s)) + dst%tau_back(s,me) if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal enddo @@ -1076,20 +1073,20 @@ module subroutine plastic_nonlocal_dotState(Mp, Temperature,timestep, & isBCC: if (lattice_structure(ph) == LATTICE_bcc_ID) then forall (s = 1:ns, sum(abs(v(s,1:4))) > 0.0_pReal) rhoDotMultiplication(s,1:2) = sum(abs(gdot(s,3:4))) / prm%b_sl(s) & ! assuming double-cross-slip of screws to be decisive for multiplication - * sqrt(stt%rho_forest(s,of)) / prm%i_sl(s) ! & ! mean free path + * sqrt(stt%rho_forest(s,me)) / prm%i_sl(s) ! & ! mean free path ! * 2.0_pReal * sum(abs(v(s,3:4))) / sum(abs(v(s,1:4))) ! ratio of screw to overall velocity determines edge generation rhoDotMultiplication(s,3:4) = sum(abs(gdot(s,3:4))) /prm%b_sl(s) & ! assuming double-cross-slip of screws to be decisive for multiplication - * sqrt(stt%rho_forest(s,of)) / prm%i_sl(s) ! & ! mean free path + * sqrt(stt%rho_forest(s,me)) / prm%i_sl(s) ! & ! mean free path ! * 2.0_pReal * sum(abs(v(s,1:2))) / sum(abs(v(s,1:4))) ! ratio of edge to overall velocity determines screw generation endforall else isBCC rhoDotMultiplication(:,1:4) = spread( & (sum(abs(gdot(:,1:2)),2) * prm%f_ed_mult + sum(abs(gdot(:,3:4)),2)) & - * sqrt(stt%rho_forest(:,of)) / prm%i_sl / prm%b_sl, 2, 4) + * sqrt(stt%rho_forest(:,me)) / prm%i_sl / prm%b_sl, 2, 4) endif isBCC - forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,instance),of) + forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,phase_plasticityInstance(ph)),me) !**************************************************************************** @@ -1132,10 +1129,10 @@ module subroutine plastic_nonlocal_dotState(Mp, Temperature,timestep, & if (lattice_structure(ph) == LATTICE_fcc_ID) & forall (s = 1:ns, prm%colinearSystem(s) > 0) & rhoDotAthermalAnnihilation(prm%colinearSystem(s),1:2) = - rhoDotAthermalAnnihilation(s,10) & - * 0.25_pReal * sqrt(stt%rho_forest(s,of)) * (dUpper(s,2) + dLower(s,2)) * prm%f_ed + * 0.25_pReal * sqrt(stt%rho_forest(s,me)) * (dUpper(s,2) + dLower(s,2)) * prm%f_ed - !*** thermally activated annihilation of edge dipoles by climb + !*** thermally activated annihilation me edge dipoles by climb rhoDotThermalAnnihilation = 0.0_pReal selfDiffusion = prm%D_0 * exp(-prm%Q_cl / (kB * Temperature)) vClimb = prm%V_at * selfDiffusion * prm%mu & @@ -1145,7 +1142,7 @@ module subroutine plastic_nonlocal_dotState(Mp, Temperature,timestep, & - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) & - rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have - rhoDot = rhoDotFlux(timestep, instance,of,ip,el) & + rhoDot = rhoDotFlux(timestep, phase_plasticityInstance(ph),me,ip,el) & + rhoDotMultiplication & + rhoDotSingle2DipoleGlide & + rhoDotAthermalAnnihilation & @@ -1162,25 +1159,25 @@ module subroutine plastic_nonlocal_dotState(Mp, Temperature,timestep, & #endif plasticState(ph)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN) else - dot%rho(:,of) = pack(rhoDot,.true.) - dot%gamma(:,of) = sum(gdot,2) + dot%rho(:,me) = pack(rhoDot,.true.) + dot%gamma(:,me) = sum(gdot,2) endif end associate -end subroutine plastic_nonlocal_dotState +end subroutine nonlocal_dotState !--------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !--------------------------------------------------------------------------------------------------- -function rhoDotFlux(timestep,instance,of,ip,el) +function rhoDotFlux(timestep,instance,me,ip,el) real(pReal), intent(in) :: & timestep !< substepped crystallite time increment integer, intent(in) :: & instance, & - of, & + me, & ip, & !< current integration point el !< current element number @@ -1243,16 +1240,16 @@ function rhoDotFlux(timestep,instance,of,ip,el) gdot = 0.0_pReal - rho = getRho(instance,of,ip,el) + rho = getRho(instance,me,ip,el) rhoSgl = rho(:,sgl) - rho0 = getRho0(instance,of,ip,el) + rho0 = getRho0(instance,me,ip,el) my_rhoSgl0 = rho0(:,sgl) - forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,instance),of) !ToDo: MD: I think we should use state0 here + forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,instance),me) !ToDo: MD: I think we should use state0 here gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4) - forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,instance),of) + forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,instance),me) !**************************************************************************** !*** calculate dislocation fluxes (only for nonlocal plasticity) @@ -1287,8 +1284,8 @@ function rhoDotFlux(timestep,instance,of,ip,el) m(1:3,:,3) = -prm%slip_transverse m(1:3,:,4) = prm%slip_transverse - my_F = constitutive_mech_F(ph)%data(1:3,1:3,of) - my_Fe = matmul(my_F, math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,of))) + my_F = constitutive_mech_F(ph)%data(1:3,1:3,me) + my_Fe = matmul(my_F, math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,me))) neighbors: do n = 1,nIPneighbors @@ -1789,14 +1786,14 @@ end subroutine kinetics !> @brief returns copy of current dislocation densities from state !> @details raw values is rectified !-------------------------------------------------------------------------------------------------- -pure function getRho(instance,of,ip,el) +pure function getRho(instance,me,ip,el) - integer, intent(in) :: instance, of,ip,el + integer, intent(in) :: instance, me,ip,el real(pReal), dimension(param(instance)%sum_N_sl,10) :: getRho associate(prm => param(instance)) - getRho = reshape(state(instance)%rho(:,of),[prm%sum_N_sl,10]) + getRho = reshape(state(instance)%rho(:,me),[prm%sum_N_sl,10]) ! ensure positive densities (not for imm, they have a sign) getRho(:,mob) = max(getRho(:,mob),0.0_pReal) @@ -1814,14 +1811,14 @@ end function getRho !> @brief returns copy of current dislocation densities from state !> @details raw values is rectified !-------------------------------------------------------------------------------------------------- -pure function getRho0(instance,of,ip,el) +pure function getRho0(instance,me,ip,el) - integer, intent(in) :: instance, of,ip,el + integer, intent(in) :: instance, me,ip,el real(pReal), dimension(param(instance)%sum_N_sl,10) :: getRho0 associate(prm => param(instance)) - getRho0 = reshape(state0(instance)%rho(:,of),[prm%sum_N_sl,10]) + getRho0 = reshape(state0(instance)%rho(:,me),[prm%sum_N_sl,10]) ! ensure positive densities (not for imm, they have a sign) getRho0(:,mob) = max(getRho0(:,mob),0.0_pReal) @@ -1834,4 +1831,4 @@ pure function getRho0(instance,of,ip,el) end function getRho0 -end submodule plastic_nonlocal +end submodule nonlocal diff --git a/src/constitutive_plastic_phenopowerlaw.f90 b/src/phase_mechanics_plastic_phenopowerlaw.f90 similarity index 91% rename from src/constitutive_plastic_phenopowerlaw.f90 rename to src/phase_mechanics_plastic_phenopowerlaw.f90 index 678acad27..ea2f53100 100644 --- a/src/constitutive_plastic_phenopowerlaw.f90 +++ b/src/phase_mechanics_plastic_phenopowerlaw.f90 @@ -4,7 +4,7 @@ !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief phenomenological crystal plasticity formulation using a powerlaw fitting !-------------------------------------------------------------------------------------------------- -submodule(constitutive:constitutive_mech) plastic_phenopowerlaw +submodule(phase:plastic) phenopowerlaw type :: tParameters real(pReal) :: & @@ -89,7 +89,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) mech, & pl - print'(/,a)', ' <<<+- plastic_phenopowerlaw init -+>>>' + print'(/,a)', ' <<<+- phase:mechanics:plastic:phenopowerlaw init -+>>>' myPlasticity = plastic_active('phenopowerlaw') Ninstances = count(myPlasticity) @@ -225,7 +225,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! allocate state arrays - Nconstituents = count(material_phaseAt == p) * discretization_nIPs + Nconstituents = count(material_phaseAt2 == p) sizeDotState = size(['xi_sl ','gamma_sl']) * prm%sum_N_sl & + size(['xi_tw ','gamma_tw']) * prm%sum_N_tw sizeState = sizeDotState @@ -285,7 +285,7 @@ end function plastic_phenopowerlaw_init !> @details asummes that deformation by dislocation glide affects twinned and untwinned volume ! equally (Taylor assumption). Twinning happens only in untwinned volume !-------------------------------------------------------------------------------------------------- -pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) +pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient @@ -295,23 +295,23 @@ pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,insta real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & - instance, & - of + ph, & + me integer :: & i,k,l,m,n - real(pReal), dimension(param(instance)%sum_N_sl) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & gdot_slip_pos,gdot_slip_neg, & dgdot_dtauslip_pos,dgdot_dtauslip_neg - real(pReal), dimension(param(instance)%sum_N_tw) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_tw) :: & gdot_twin,dgdot_dtautwin Lp = 0.0_pReal dLp_dMp = 0.0_pReal - associate(prm => param(instance)) + associate(prm => param(phase_plasticityInstance(ph))) - call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg) + call kinetics_slip(Mp,phase_plasticityInstance(ph),me,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg) slipSystems: do i = 1, prm%sum_N_sl Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%P_sl(1:3,1:3,i) forall (k=1:3,l=1:3,m=1:3,n=1:3) & @@ -320,7 +320,7 @@ pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,insta + dgdot_dtauslip_neg(i) * prm%P_sl(k,l,i) * prm%nonSchmid_neg(m,n,i) enddo slipSystems - call kinetics_twin(Mp,instance,of,gdot_twin,dgdot_dtautwin) + call kinetics_twin(Mp,phase_plasticityInstance(ph),me,gdot_twin,dgdot_dtautwin) twinSystems: do i = 1, prm%sum_N_tw Lp = Lp + gdot_twin(i)*prm%P_tw(1:3,1:3,i) forall (k=1:3,l=1:3,m=1:3,n=1:3) & @@ -330,32 +330,33 @@ pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,insta end associate -end subroutine plastic_phenopowerlaw_LpAndItsTangent +end subroutine phenopowerlaw_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- -module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) +module subroutine phenopowerlaw_dotState(Mp,ph,me) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & - instance, & - of + ph, & + me real(pReal) :: & c_SlipSlip,c_TwinSlip,c_TwinTwin, & xi_slip_sat_offset,& sumGamma,sumF - real(pReal), dimension(param(instance)%sum_N_sl) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & left_SlipSlip,right_SlipSlip, & gdot_slip_pos,gdot_slip_neg - associate(prm => param(instance), stt => state(instance), dot => dotState(instance)) + associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph)), & + dot => dotState(phase_plasticityInstance(ph))) - sumGamma = sum(stt%gamma_slip(:,of)) - sumF = sum(stt%gamma_twin(:,of)/prm%gamma_tw_char) + sumGamma = sum(stt%gamma_slip(:,me)) + sumF = sum(stt%gamma_twin(:,me)/prm%gamma_tw_char) !-------------------------------------------------------------------------------------------------- ! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices @@ -367,26 +368,26 @@ module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) ! calculate left and right vectors left_SlipSlip = 1.0_pReal + prm%h_int xi_slip_sat_offset = prm%f_sl_sat_tw*sqrt(sumF) - right_SlipSlip = abs(1.0_pReal-stt%xi_slip(:,of) / (prm%xi_inf_sl+xi_slip_sat_offset)) **prm%a_sl & - * sign(1.0_pReal,1.0_pReal-stt%xi_slip(:,of) / (prm%xi_inf_sl+xi_slip_sat_offset)) + right_SlipSlip = abs(1.0_pReal-stt%xi_slip(:,me) / (prm%xi_inf_sl+xi_slip_sat_offset)) **prm%a_sl & + * sign(1.0_pReal,1.0_pReal-stt%xi_slip(:,me) / (prm%xi_inf_sl+xi_slip_sat_offset)) !-------------------------------------------------------------------------------------------------- ! shear rates - call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg) - dot%gamma_slip(:,of) = abs(gdot_slip_pos+gdot_slip_neg) - call kinetics_twin(Mp,instance,of,dot%gamma_twin(:,of)) + call kinetics_slip(Mp,phase_plasticityInstance(ph),me,gdot_slip_pos,gdot_slip_neg) + dot%gamma_slip(:,me) = abs(gdot_slip_pos+gdot_slip_neg) + call kinetics_twin(Mp,phase_plasticityInstance(ph),me,dot%gamma_twin(:,me)) !-------------------------------------------------------------------------------------------------- ! hardening - dot%xi_slip(:,of) = c_SlipSlip * left_SlipSlip * & - matmul(prm%h_sl_sl,dot%gamma_slip(:,of)*right_SlipSlip) & - + matmul(prm%h_sl_tw,dot%gamma_twin(:,of)) + dot%xi_slip(:,me) = c_SlipSlip * left_SlipSlip * & + matmul(prm%h_sl_sl,dot%gamma_slip(:,me)*right_SlipSlip) & + + matmul(prm%h_sl_tw,dot%gamma_twin(:,me)) - dot%xi_twin(:,of) = c_TwinSlip * matmul(prm%h_tw_sl,dot%gamma_slip(:,of)) & - + c_TwinTwin * matmul(prm%h_tw_tw,dot%gamma_twin(:,of)) + dot%xi_twin(:,me) = c_TwinSlip * matmul(prm%h_tw_sl,dot%gamma_slip(:,me)) & + + c_TwinTwin * matmul(prm%h_tw_tw,dot%gamma_twin(:,me)) end associate -end subroutine plastic_phenopowerlaw_dotState +end subroutine phenopowerlaw_dotState !-------------------------------------------------------------------------------------------------- @@ -431,14 +432,14 @@ end subroutine plastic_phenopowerlaw_results ! 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(Mp,instance,of, & +pure subroutine kinetics_slip(Mp,instance,me, & gdot_slip_pos,gdot_slip_neg,dgdot_dtau_slip_pos,dgdot_dtau_slip_neg) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & instance, & - of + me real(pReal), intent(out), dimension(param(instance)%sum_N_sl) :: & gdot_slip_pos, & @@ -462,14 +463,14 @@ pure subroutine kinetics_slip(Mp,instance,of, & where(dNeq0(tau_slip_pos)) gdot_slip_pos = prm%dot_gamma_0_sl * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active - * sign(abs(tau_slip_pos/stt%xi_slip(:,of))**prm%n_sl, tau_slip_pos) + * sign(abs(tau_slip_pos/stt%xi_slip(:,me))**prm%n_sl, tau_slip_pos) else where gdot_slip_pos = 0.0_pReal end where where(dNeq0(tau_slip_neg)) gdot_slip_neg = prm%dot_gamma_0_sl * 0.5_pReal & ! only used if non-Schmid active, always 1/2 - * sign(abs(tau_slip_neg/stt%xi_slip(:,of))**prm%n_sl, tau_slip_neg) + * sign(abs(tau_slip_neg/stt%xi_slip(:,me))**prm%n_sl, tau_slip_neg) else where gdot_slip_neg = 0.0_pReal end where @@ -500,14 +501,14 @@ end subroutine kinetics_slip ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- -pure subroutine kinetics_twin(Mp,instance,of,& +pure subroutine kinetics_twin(Mp,instance,me,& gdot_twin,dgdot_dtau_twin) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & instance, & - of + me real(pReal), dimension(param(instance)%sum_N_tw), intent(out) :: & gdot_twin @@ -525,8 +526,8 @@ pure subroutine kinetics_twin(Mp,instance,of,& enddo where(tau_twin > 0.0_pReal) - gdot_twin = (1.0_pReal-sum(stt%gamma_twin(:,of)/prm%gamma_tw_char)) & ! only twin in untwinned volume fraction - * prm%dot_gamma_0_tw*(abs(tau_twin)/stt%xi_twin(:,of))**prm%n_tw + gdot_twin = (1.0_pReal-sum(stt%gamma_twin(:,me)/prm%gamma_tw_char)) & ! only twin in untwinned volume fraction + * prm%dot_gamma_0_tw*(abs(tau_twin)/stt%xi_twin(:,me))**prm%n_tw else where gdot_twin = 0.0_pReal end where @@ -543,4 +544,4 @@ pure subroutine kinetics_twin(Mp,instance,of,& end subroutine kinetics_twin -end submodule plastic_phenopowerlaw +end submodule phenopowerlaw diff --git a/src/constitutive_thermal.f90 b/src/phase_thermal.f90 similarity index 67% rename from src/constitutive_thermal.f90 rename to src/phase_thermal.f90 index 831fd236d..d5d52b010 100644 --- a/src/constitutive_thermal.f90 +++ b/src/phase_thermal.f90 @@ -1,7 +1,7 @@ !---------------------------------------------------------------------------------------------------- !> @brief internal microstructure state for all thermal sources and kinematics constitutive models !---------------------------------------------------------------------------------------------------- -submodule(constitutive) constitutive_thermal +submodule(phase) thermal enum, bind(c); enumerator :: & THERMAL_UNDEFINED_ID ,& @@ -18,49 +18,44 @@ submodule(constitutive) constitutive_thermal type(tDataContainer), dimension(:), allocatable :: current integer :: thermal_source_maxSizeDotState + + interface - module function source_thermal_dissipation_init(source_length) result(mySources) - integer, intent(in) :: source_length - logical, dimension(:,:), allocatable :: mySources - end function source_thermal_dissipation_init + module function dissipation_init(source_length) result(mySources) + integer, intent(in) :: source_length + logical, dimension(:,:), allocatable :: mySources + end function dissipation_init - module function source_thermal_externalheat_init(source_length) result(mySources) - integer, intent(in) :: source_length - logical, dimension(:,:), allocatable :: mySources - end function source_thermal_externalheat_init - - module function kinematics_thermal_expansion_init(kinematics_length) result(myKinematics) - integer, intent(in) :: kinematics_length - logical, dimension(:,:), allocatable :: myKinematics - end function kinematics_thermal_expansion_init + module function externalheat_init(source_length) result(mySources) + integer, intent(in) :: source_length + logical, dimension(:,:), allocatable :: mySources + end function externalheat_init - module subroutine source_thermal_externalheat_dotState(ph, me) + + + module subroutine externalheat_dotState(ph, me) integer, intent(in) :: & ph, & me - end subroutine source_thermal_externalheat_dotState + end subroutine externalheat_dotState + module subroutine dissipation_getRate(TDot, ph,me) + integer, intent(in) :: & + ph, & + me + real(pReal), intent(out) :: & + TDot + end subroutine dissipation_getRate - module subroutine thermal_dissipation_getRate(TDot, Tstar,Lp,phase) - integer, intent(in) :: & - phase !< phase ID of element - real(pReal), intent(in), dimension(3,3) :: & - Tstar !< 2nd Piola Kirchhoff stress tensor for a given element - real(pReal), intent(in), dimension(3,3) :: & - Lp !< plastic velocuty gradient for a given element - real(pReal), intent(out) :: & - TDot - end subroutine thermal_dissipation_getRate - - module subroutine thermal_externalheat_getRate(TDot, ph,me) - integer, intent(in) :: & - ph, & - me - real(pReal), intent(out) :: & - TDot - end subroutine thermal_externalheat_getRate + module subroutine externalheat_getRate(TDot, ph,me) + integer, intent(in) :: & + ph, & + me + real(pReal), intent(out) :: & + TDot + end subroutine externalheat_getRate end interface @@ -82,7 +77,7 @@ module subroutine thermal_init(phases) Nconstituents - print'(/,a)', ' <<<+- constitutive_thermal init -+>>>' + print'(/,a)', ' <<<+- phase:thermal init -+>>>' allocate(current(phases%length)) @@ -91,7 +86,7 @@ module subroutine thermal_init(phases) do ph = 1, phases%length - Nconstituents = count(material_phaseAt == ph) * discretization_nIPs + Nconstituents = count(material_phaseAt2 == ph) allocate(current(ph)%T(Nconstituents),source=300.0_pReal) allocate(current(ph)%dot_T(Nconstituents),source=0.0_pReal) @@ -108,8 +103,8 @@ module subroutine thermal_init(phases) allocate(thermal_source(maxval(thermal_Nsources),phases%length), source = THERMAL_UNDEFINED_ID) if(maxval(thermal_Nsources) /= 0) then - where(source_thermal_dissipation_init (maxval(thermal_Nsources))) thermal_source = THERMAL_DISSIPATION_ID - where(source_thermal_externalheat_init(maxval(thermal_Nsources))) thermal_source = THERMAL_EXTERNALHEAT_ID + where(dissipation_init (maxval(thermal_Nsources))) thermal_source = THERMAL_DISSIPATION_ID + where(externalheat_init(maxval(thermal_Nsources))) thermal_source = THERMAL_EXTERNALHEAT_ID endif thermal_source_maxSizeDotState = 0 @@ -123,54 +118,40 @@ module subroutine thermal_init(phases) maxval(thermalState(ph)%p%sizeDotState)) enddo PhaseLoop2 -!-------------------------------------------------------------------------------------------------- -!initialize kinematic mechanisms - if(maxval(phase_Nkinematics) /= 0) where(kinematics_thermal_expansion_init(maxval(phase_Nkinematics))) & - phase_kinematics = KINEMATICS_thermal_expansion_ID - end subroutine thermal_init !---------------------------------------------------------------------------------------------- !< @brief calculates thermal dissipation rate !---------------------------------------------------------------------------------------------- -module subroutine constitutive_thermal_getRate(TDot, ip, el) +module subroutine constitutive_thermal_getRate(TDot, ph,me) - integer, intent(in) :: & - ip, & !< integration point number - el !< element number + integer, intent(in) :: ph, me real(pReal), intent(out) :: & TDot real(pReal) :: & my_Tdot integer :: & - ph, & - homog, & - me, & - so, & - co + so - homog = material_homogenizationAt(el) TDot = 0.0_pReal - do co = 1, homogenization_Nconstituents(homog) - ph = material_phaseAt(co,el) - me = material_phasememberAt(co,ip,el) - do so = 1, thermal_Nsources(ph) - select case(thermal_source(so,ph)) - case (THERMAL_DISSIPATION_ID) - call thermal_dissipation_getRate(my_Tdot, mech_S(ph,me),mech_L_p(ph,me),ph) - case (THERMAL_EXTERNALHEAT_ID) - call thermal_externalheat_getRate(my_Tdot, ph,me) + do so = 1, thermal_Nsources(ph) + select case(thermal_source(so,ph)) + case (THERMAL_DISSIPATION_ID) + call dissipation_getRate(my_Tdot, ph,me) + + case (THERMAL_EXTERNALHEAT_ID) + call externalheat_getRate(my_Tdot, ph,me) + + case default + my_Tdot = 0.0_pReal + end select + Tdot = Tdot + my_Tdot + enddo - case default - my_Tdot = 0.0_pReal - end select - Tdot = Tdot + my_Tdot - enddo - enddo end subroutine constitutive_thermal_getRate @@ -191,7 +172,7 @@ function constitutive_thermal_collectDotState(ph,me) result(broken) SourceLoop: do i = 1, thermal_Nsources(ph) if (thermal_source(i,ph) == THERMAL_EXTERNALHEAT_ID) & - call source_thermal_externalheat_dotState(ph,me) + call externalheat_dotState(ph,me) broken = broken .or. any(IEEE_is_NaN(thermalState(ph)%p(i)%dotState(:,me))) @@ -206,8 +187,6 @@ module function thermal_stress(Delta_t,ph,me) result(converged_) integer, intent(in) :: ph, me logical :: converged_ - integer :: so - converged_ = .not. integrateThermalState(Delta_t,ph,me) @@ -268,6 +247,20 @@ module function thermal_T(ph,me) result(T) end function thermal_T +!---------------------------------------------------------------------------------------------- +!< @brief Get rate of temperature (for use by non-thermal physics) +!---------------------------------------------------------------------------------------------- +module function thermal_dot_T(ph,me) result(dot_T) + + integer, intent(in) :: ph, me + real(pReal) :: dot_T + + + dot_T = current(ph)%dot_T(me) + +end function thermal_dot_T + + !---------------------------------------------------------------------------------------------- !< @brief Set temperature !---------------------------------------------------------------------------------------------- @@ -318,4 +311,4 @@ function thermal_active(source_label,src_length) result(active_source) end function thermal_active -end submodule constitutive_thermal +end submodule thermal diff --git a/src/constitutive_thermal_dissipation.f90 b/src/phase_thermal_dissipation.f90 similarity index 52% rename from src/constitutive_thermal_dissipation.f90 rename to src/phase_thermal_dissipation.f90 index ae2d5735e..7a043308d 100644 --- a/src/constitutive_thermal_dissipation.f90 +++ b/src/phase_thermal_dissipation.f90 @@ -4,11 +4,7 @@ !> @brief material subroutine for thermal source due to plastic dissipation !> @details to be done !-------------------------------------------------------------------------------------------------- -submodule(constitutive:constitutive_thermal) source_dissipation - - integer, dimension(:), allocatable :: & - source_thermal_dissipation_offset, & !< which source is my current thermal dissipation mechanism? - source_thermal_dissipation_instance !< instance of thermal dissipation source mechanism +submodule(phase:thermal) dissipation type :: tParameters !< container type for internal constitutive parameters real(pReal) :: & @@ -25,7 +21,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -module function source_thermal_dissipation_init(source_length) result(mySources) +module function dissipation_init(source_length) result(mySources) integer, intent(in) :: source_length logical, dimension(:,:), allocatable :: mySources @@ -35,9 +31,9 @@ module function source_thermal_dissipation_init(source_length) result(mySources) phase, & sources, thermal, & src - integer :: Ninstances,sourceOffset,Nconstituents,p + integer :: Ninstances,so,Nconstituents,ph - print'(/,a)', ' <<<+- thermal_dissipation init -+>>>' + print'(/,a)', ' <<<+- phase:thermal:dissipation init -+>>>' mySources = thermal_active('dissipation',source_length) @@ -46,25 +42,21 @@ module function source_thermal_dissipation_init(source_length) result(mySources) if(Ninstances == 0) return phases => config_material%get('phase') - allocate(param(Ninstances)) - allocate(source_thermal_dissipation_offset (phases%length), source=0) - allocate(source_thermal_dissipation_instance(phases%length), source=0) + allocate(param(phases%length)) - do p = 1, phases%length - phase => phases%get(p) - if(any(mySources(:,p))) source_thermal_dissipation_instance(p) = count(mySources(:,1:p)) - if(count(mySources(:,p)) == 0) cycle + do ph = 1, phases%length + phase => phases%get(ph) + if(count(mySources(:,ph)) == 0) cycle !ToDo: error if > 1 thermal => phase%get('thermal') sources => thermal%get('source') - do sourceOffset = 1, sources%length - if(mySources(sourceOffset,p)) then - source_thermal_dissipation_offset(p) = sourceOffset - associate(prm => param(source_thermal_dissipation_instance(p))) - src => sources%get(sourceOffset) + do so = 1, sources%length + if(mySources(so,ph)) then + associate(prm => param(ph)) + src => sources%get(so) - prm%kappa = src%get_asFloat('kappa') - Nconstituents = count(material_phaseAt==p) * discretization_nIPs - call constitutive_allocateState(thermalState(p)%p(sourceOffset),Nconstituents,0,0,0) + prm%kappa = src%get_asFloat('kappa') + Nconstituents = count(material_phaseAt2 == ph) + call constitutive_allocateState(thermalState(ph)%p(so),Nconstituents,0,0,0) end associate endif @@ -72,28 +64,23 @@ module function source_thermal_dissipation_init(source_length) result(mySources) enddo -end function source_thermal_dissipation_init +end function dissipation_init !-------------------------------------------------------------------------------------------------- !> @brief Ninstancess dissipation rate !-------------------------------------------------------------------------------------------------- -module subroutine thermal_dissipation_getRate(TDot, Tstar, Lp, phase) - - integer, intent(in) :: & - phase - real(pReal), intent(in), dimension(3,3) :: & - Tstar - real(pReal), intent(in), dimension(3,3) :: & - Lp +module subroutine dissipation_getRate(TDot, ph,me) + integer, intent(in) :: ph, me real(pReal), intent(out) :: & TDot - associate(prm => param(source_thermal_dissipation_instance(phase))) - TDot = prm%kappa*sum(abs(Tstar*Lp)) + + associate(prm => param(ph)) + TDot = prm%kappa*sum(abs(mech_S(ph,me)*mech_L_p(ph,me))) end associate -end subroutine thermal_dissipation_getRate +end subroutine dissipation_getRate -end submodule source_dissipation +end submodule dissipation diff --git a/src/phase_thermal_externalheat.f90 b/src/phase_thermal_externalheat.f90 new file mode 100644 index 000000000..ddae6763b --- /dev/null +++ b/src/phase_thermal_externalheat.f90 @@ -0,0 +1,134 @@ +!-------------------------------------------------------------------------------------------------- +!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH +!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH +!> @author Philip Eisenlohr, Michigan State University +!> @brief material subroutine for variable heat source +!-------------------------------------------------------------------------------------------------- +submodule(phase:thermal) externalheat + + + integer, dimension(:), allocatable :: & + source_thermal_externalheat_offset !< which source is my current thermal dissipation mechanism? + + type :: tParameters !< container type for internal constitutive parameters + real(pReal), dimension(:), allocatable :: & + t_n, & + f_T + integer :: & + nIntervals + end type tParameters + + type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances) + + +contains + + +!-------------------------------------------------------------------------------------------------- +!> @brief module initialization +!> @details reads in material parameters, allocates arrays, and does sanity checks +!-------------------------------------------------------------------------------------------------- +module function externalheat_init(source_length) result(mySources) + + integer, intent(in) :: source_length + logical, dimension(:,:), allocatable :: mySources + + class(tNode), pointer :: & + phases, & + phase, & + sources, thermal, & + src + integer :: Ninstances,so,Nconstituents,ph + + print'(/,a)', ' <<<+- phase:thermal:externalheat init -+>>>' + + mySources = thermal_active('externalheat',source_length) + + Ninstances = count(mySources) + print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) + if(Ninstances == 0) return + + phases => config_material%get('phase') + allocate(param(phases%length)) + allocate(source_thermal_externalheat_offset (phases%length), source=0) + + do ph = 1, phases%length + phase => phases%get(ph) + if(count(mySources(:,ph)) == 0) cycle + thermal => phase%get('thermal') + sources => thermal%get('source') + do so = 1, sources%length + if(mySources(so,ph)) then + source_thermal_externalheat_offset(ph) = so + associate(prm => param(ph)) + src => sources%get(so) + + prm%t_n = src%get_asFloats('t_n') + prm%nIntervals = size(prm%t_n) - 1 + + prm%f_T = src%get_asFloats('f_T',requiredSize = size(prm%t_n)) + + Nconstituents = count(material_phaseAt2 == ph) + call constitutive_allocateState(thermalState(ph)%p(so),Nconstituents,1,1,0) + end associate + endif + enddo + enddo + +end function externalheat_init + + +!-------------------------------------------------------------------------------------------------- +!> @brief rate of change of state +!> @details state only contains current time to linearly interpolate given heat powers +!-------------------------------------------------------------------------------------------------- +module subroutine externalheat_dotState(ph, me) + + integer, intent(in) :: & + ph, & + me + + integer :: & + so + + so = source_thermal_externalheat_offset(ph) + + thermalState(ph)%p(so)%dotState(1,me) = 1.0_pReal ! state is current time + +end subroutine externalheat_dotState + + +!-------------------------------------------------------------------------------------------------- +!> @brief returns local heat generation rate +!-------------------------------------------------------------------------------------------------- +module subroutine externalheat_getRate(TDot, ph, me) + + integer, intent(in) :: & + ph, & + me + real(pReal), intent(out) :: & + TDot + + integer :: & + so, interval + real(pReal) :: & + frac_time + + so = source_thermal_externalheat_offset(ph) + + associate(prm => param(ph)) + do interval = 1, prm%nIntervals ! scan through all rate segments + frac_time = (thermalState(ph)%p(so)%state(1,me) - prm%t_n(interval)) & + / (prm%t_n(interval+1) - prm%t_n(interval)) ! fractional time within segment + if ( (frac_time < 0.0_pReal .and. interval == 1) & + .or. (frac_time >= 1.0_pReal .and. interval == prm%nIntervals) & + .or. (frac_time >= 0.0_pReal .and. frac_time < 1.0_pReal) ) & + TDot = prm%f_T(interval ) * (1.0_pReal - frac_time) + & + prm%f_T(interval+1) * frac_time ! interpolate heat rate between segment boundaries... + ! ...or extrapolate if outside me bounds + enddo + end associate + +end subroutine externalheat_getRate + +end submodule externalheat