From e52a7477742a8f59dd0e5e9ab3051aa2e11dc22d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 2 Dec 2019 20:43:02 +0100 Subject: [PATCH 01/20] submodules allow inter-module communication --- src/constitutive.f90 | 87 +++++++++++++++++++++++++++++++++-- src/plastic_isotropic.f90 | 32 +++---------- src/plastic_none.f90 | 15 ++---- src/plastic_phenopowerlaw.f90 | 30 +++--------- 4 files changed, 100 insertions(+), 64 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 4067c026a..54dbe39fa 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -4,6 +4,7 @@ !> @brief elasticity, plasticity, internal microstructure state !-------------------------------------------------------------------------------------------------- module constitutive + use prec use math use debug use numerics @@ -14,9 +15,6 @@ module constitutive use HDF5_utilities use lattice use discretization - use plastic_none - use plastic_isotropic - use plastic_phenopowerlaw use plastic_kinehardening use plastic_dislotwin use plastic_disloucla @@ -40,6 +38,89 @@ module constitutive constitutive_source_maxSizePostResults, & constitutive_source_maxSizeDotState + + interface + + module subroutine plastic_none_init + end subroutine plastic_none_init + + module subroutine plastic_isotropic_init + end subroutine plastic_isotropic_init + + module subroutine plastic_phenopowerlaw_init + end subroutine plastic_phenopowerlaw_init + + + module subroutine plastic_isotropic_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_isotropic_LpAndItsTangent + + 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 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 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_isotropic_results(instance,group) + integer, intent(in) :: instance + character(len=*), intent(in) :: group + end subroutine plastic_isotropic_results + + module subroutine plastic_phenopowerlaw_results(instance,group) + integer, intent(in) :: instance + character(len=*), intent(in) :: group + end subroutine plastic_phenopowerlaw_results + + end interface + public :: & constitutive_init, & constitutive_homogenizedC, & diff --git a/src/plastic_isotropic.f90 b/src/plastic_isotropic.f90 index 9beb2262b..86b764022 100644 --- a/src/plastic_isotropic.f90 +++ b/src/plastic_isotropic.f90 @@ -7,18 +7,7 @@ !! resolving the stress on the slip systems. Will give the response of phenopowerlaw for an !! untextured polycrystal !-------------------------------------------------------------------------------------------------- -module plastic_isotropic - use prec - use debug - use math - use IO - use material - use config - use discretization - use results - - implicit none - private +submodule(constitutive) plastic_isotropic enum, bind(c) enumerator :: & @@ -64,20 +53,13 @@ module plastic_isotropic dotState, & state - public :: & - plastic_isotropic_init, & - plastic_isotropic_LpAndItsTangent, & - plastic_isotropic_LiAndItsTangent, & - plastic_isotropic_dotState, & - plastic_isotropic_results - contains !-------------------------------------------------------------------------------------------------- !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine plastic_isotropic_init +module subroutine plastic_isotropic_init integer :: & Ninstance, & @@ -210,7 +192,7 @@ end subroutine plastic_isotropic_init !-------------------------------------------------------------------------------------------------- !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- -subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) +module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient @@ -271,7 +253,7 @@ end subroutine plastic_isotropic_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- -subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of) +module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of) real(pReal), dimension(3,3), intent(out) :: & Li !< inleastic velocity gradient @@ -323,7 +305,7 @@ subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of) !-------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -subroutine plastic_isotropic_dotState(Mp,instance,of) +module subroutine plastic_isotropic_dotState(Mp,instance,of) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -372,7 +354,7 @@ end subroutine plastic_isotropic_dotState !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- -subroutine plastic_isotropic_results(instance,group) +module subroutine plastic_isotropic_results(instance,group) #if defined(PETSc) || defined(DAMASK_HDF5) integer, intent(in) :: instance @@ -396,4 +378,4 @@ subroutine plastic_isotropic_results(instance,group) end subroutine plastic_isotropic_results -end module plastic_isotropic +end submodule plastic_isotropic diff --git a/src/plastic_none.f90 b/src/plastic_none.f90 index a4979bb2c..14f728cc1 100644 --- a/src/plastic_none.f90 +++ b/src/plastic_none.f90 @@ -4,16 +4,7 @@ !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief Dummy plasticity for purely elastic material !-------------------------------------------------------------------------------------------------- -module plastic_none - use material - use discretization - use debug - - implicit none - private - - public :: & - plastic_none_init +submodule(constitutive) plastic_none contains @@ -21,7 +12,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine plastic_none_init +module subroutine plastic_none_init integer :: & Ninstance, & @@ -44,4 +35,4 @@ subroutine plastic_none_init end subroutine plastic_none_init -end module plastic_none +end submodule plastic_none diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 index b8f5c8306..aa2f47435 100644 --- a/src/plastic_phenopowerlaw.f90 +++ b/src/plastic_phenopowerlaw.f90 @@ -4,19 +4,7 @@ !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief phenomenological crystal plasticity formulation using a powerlaw fitting !-------------------------------------------------------------------------------------------------- -module plastic_phenopowerlaw - use prec - use debug - use math - use IO - use material - use config - use lattice - use discretization - use results - - implicit none - private +submodule(constitutive) plastic_phenopowerlaw enum, bind(c) enumerator :: & @@ -91,12 +79,6 @@ module plastic_phenopowerlaw dotState, & state - public :: & - plastic_phenopowerlaw_init, & - plastic_phenopowerlaw_LpAndItsTangent, & - plastic_phenopowerlaw_dotState, & - plastic_phenopowerlaw_results - contains @@ -104,7 +86,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine plastic_phenopowerlaw_init +module subroutine plastic_phenopowerlaw_init integer :: & Ninstance, & @@ -356,7 +338,7 @@ end subroutine 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 subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) +pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient @@ -407,7 +389,7 @@ end subroutine plastic_phenopowerlaw_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) +module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -463,7 +445,7 @@ end subroutine plastic_phenopowerlaw_dotState !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- -subroutine plastic_phenopowerlaw_results(instance,group) +module subroutine plastic_phenopowerlaw_results(instance,group) #if defined(PETSc) || defined(DAMASK_HDF5) integer, intent(in) :: instance @@ -622,4 +604,4 @@ pure subroutine kinetics_twin(Mp,instance,of,& end subroutine kinetics_twin -end module plastic_phenopowerlaw +end submodule plastic_phenopowerlaw From 226b715c460eb2426178950af5a2d04aca6b55da Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 2 Dec 2019 21:38:41 +0100 Subject: [PATCH 02/20] file names reflect hierarchical structure --- src/constitutive.f90 | 44 +++++++++++++++++-- ...f90 => constitutive_plastic_isotropic.f90} | 0 ...=> constitutive_plastic_kinehardening.f90} | 37 ++++------------ ...none.f90 => constitutive_plastic_none.f90} | 0 ...=> constitutive_plastic_phenopowerlaw.f90} | 0 src/plastic_dislotwin.f90 | 2 +- 6 files changed, 50 insertions(+), 33 deletions(-) rename src/{plastic_isotropic.f90 => constitutive_plastic_isotropic.f90} (100%) rename src/{plastic_kinematichardening.f90 => constitutive_plastic_kinehardening.f90} (96%) rename src/{plastic_none.f90 => constitutive_plastic_none.f90} (100%) rename src/{plastic_phenopowerlaw.f90 => constitutive_plastic_phenopowerlaw.f90} (100%) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 54dbe39fa..2549237d6 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -15,7 +15,6 @@ module constitutive use HDF5_utilities use lattice use discretization - use plastic_kinehardening use plastic_dislotwin use plastic_disloucla use plastic_nonlocal @@ -49,7 +48,9 @@ module constitutive module subroutine plastic_phenopowerlaw_init end subroutine plastic_phenopowerlaw_init - + + module subroutine plastic_kinehardening_init + end subroutine plastic_kinehardening_init module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) real(pReal), dimension(3,3), intent(out) :: & @@ -76,6 +77,19 @@ module constitutive instance, & of end subroutine plastic_phenopowerlaw_LpAndItsTangent + + 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_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of) @@ -108,6 +122,23 @@ module constitutive 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_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 plastic_isotropic_results(instance,group) integer, intent(in) :: instance @@ -117,8 +148,13 @@ module constitutive module subroutine plastic_phenopowerlaw_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group - end subroutine plastic_phenopowerlaw_results - + end subroutine plastic_phenopowerlaw_results + + module subroutine plastic_kinehardening_results(instance,group) + integer, intent(in) :: instance + character(len=*), intent(in) :: group + end subroutine plastic_kinehardening_results + end interface public :: & diff --git a/src/plastic_isotropic.f90 b/src/constitutive_plastic_isotropic.f90 similarity index 100% rename from src/plastic_isotropic.f90 rename to src/constitutive_plastic_isotropic.f90 diff --git a/src/plastic_kinematichardening.f90 b/src/constitutive_plastic_kinehardening.f90 similarity index 96% rename from src/plastic_kinematichardening.f90 rename to src/constitutive_plastic_kinehardening.f90 index 2a3dc4640..39c3fc01a 100644 --- a/src/plastic_kinematichardening.f90 +++ b/src/constitutive_plastic_kinehardening.f90 @@ -5,19 +5,7 @@ !> @brief Phenomenological crystal plasticity using a power law formulation for the shear rates !! and a Voce-type kinematic hardening rule !-------------------------------------------------------------------------------------------------- -module plastic_kinehardening - use prec - use debug - use math - use IO - use material - use config - use lattice - use discretization - use results - - implicit none - private +submodule(constitutive) plastic_kinehardening enum, bind(c) enumerator :: & @@ -80,13 +68,6 @@ module plastic_kinehardening deltaState, & state - public :: & - plastic_kinehardening_init, & - plastic_kinehardening_LpAndItsTangent, & - plastic_kinehardening_dotState, & - plastic_kinehardening_deltaState, & - plastic_kinehardening_results - contains @@ -94,7 +75,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine plastic_kinehardening_init +module subroutine plastic_kinehardening_init integer :: & Ninstance, & @@ -304,7 +285,7 @@ end subroutine plastic_kinehardening_init !-------------------------------------------------------------------------------------------------- !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- -pure subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) +pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient @@ -346,7 +327,7 @@ end subroutine plastic_kinehardening_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -subroutine plastic_kinehardening_dotState(Mp,instance,of) +module subroutine plastic_kinehardening_dotState(Mp,instance,of) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -388,7 +369,7 @@ end subroutine plastic_kinehardening_dotState !-------------------------------------------------------------------------------------------------- !> @brief calculates (instantaneous) incremental change of microstructure !-------------------------------------------------------------------------------------------------- -subroutine plastic_kinehardening_deltaState(Mp,instance,of) +module subroutine plastic_kinehardening_deltaState(Mp,instance,of) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -436,11 +417,11 @@ end subroutine plastic_kinehardening_deltaState !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- -subroutine plastic_kinehardening_results(instance,group) +module subroutine plastic_kinehardening_results(instance,group) #if defined(PETSc) || defined(DAMASK_HDF5) - integer, intent(in) :: instance - character(len=*) :: group + integer, intent(in) :: instance + character(len=*), intent(in) :: group integer :: o associate(prm => param(instance), stt => state(instance)) @@ -548,4 +529,4 @@ pure subroutine kinetics(Mp,instance,of, & end subroutine kinetics -end module plastic_kinehardening +end submodule plastic_kinehardening diff --git a/src/plastic_none.f90 b/src/constitutive_plastic_none.f90 similarity index 100% rename from src/plastic_none.f90 rename to src/constitutive_plastic_none.f90 diff --git a/src/plastic_phenopowerlaw.f90 b/src/constitutive_plastic_phenopowerlaw.f90 similarity index 100% rename from src/plastic_phenopowerlaw.f90 rename to src/constitutive_plastic_phenopowerlaw.f90 diff --git a/src/plastic_dislotwin.f90 b/src/plastic_dislotwin.f90 index 1f731a891..24b6bd3da 100644 --- a/src/plastic_dislotwin.f90 +++ b/src/plastic_dislotwin.f90 @@ -21,7 +21,7 @@ module plastic_dislotwin implicit none private - real(pReal), parameter :: & + real(pReal), parameter :: & kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin enum, bind(c) From 4ee2e551b846c5d2dc4a428fd7b268458ecaa735 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 2 Dec 2019 21:51:25 +0100 Subject: [PATCH 03/20] migrating to submodules --- src/constitutive.f90 | 100 ++++++++++++++++-- ...f90 => constitutive_plastic_disloUCLA.f90} | 59 ++++------- ...f90 => constitutive_plastic_dislotwin.f90} | 46 +++----- 3 files changed, 126 insertions(+), 79 deletions(-) rename src/{plastic_disloUCLA.f90 => constitutive_plastic_disloUCLA.f90} (95%) rename src/{plastic_dislotwin.f90 => constitutive_plastic_dislotwin.f90} (98%) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 2549237d6..7f87d3cda 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -12,11 +12,8 @@ module constitutive use config use material use results - use HDF5_utilities use lattice use discretization - use plastic_dislotwin - use plastic_disloucla use plastic_nonlocal use geometry_plastic_nonlocal use source_thermal_dissipation @@ -37,7 +34,6 @@ module constitutive constitutive_source_maxSizePostResults, & constitutive_source_maxSizeDotState - interface module subroutine plastic_none_init @@ -51,7 +47,14 @@ module constitutive module subroutine plastic_kinehardening_init end subroutine plastic_kinehardening_init + + module subroutine plastic_dislotwin_init + end subroutine plastic_dislotwin_init + module subroutine plastic_disloUCLA_init + end subroutine plastic_disloUCLA_init + + module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient @@ -91,6 +94,36 @@ module constitutive 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_disloUCLA_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_disloUCLA_LpAndItsTangent + module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of) real(pReal), dimension(3,3), intent(out) :: & @@ -129,6 +162,41 @@ module constitutive 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_disloUCLA_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_disloUCLA_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_disloUCLA_dependentState(instance,of) + integer, intent(in) :: & + instance, & + of + end subroutine plastic_disloUCLA_dependentState module subroutine plastic_kinehardening_deltaState(Mp,instance,of) @@ -138,6 +206,16 @@ module constitutive instance, & of end subroutine plastic_kinehardening_deltaState + + + module function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC) + real(pReal), dimension(6,6) :: & + homogenizedC + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + end function plastic_dislotwin_homogenizedC module subroutine plastic_isotropic_results(instance,group) @@ -154,6 +232,16 @@ module constitutive integer, intent(in) :: instance character(len=*), intent(in) :: group end subroutine plastic_kinehardening_results + + module subroutine plastic_dislotwin_results(instance,group) + integer, intent(in) :: instance + character(len=*), intent(in) :: group + end subroutine plastic_dislotwin_results + + module subroutine plastic_disloUCLA_results(instance,group) + integer, intent(in) :: instance + character(len=*), intent(in) :: group + end subroutine plastic_disloUCLA_results end interface @@ -810,11 +898,11 @@ subroutine constitutive_results character(len=256) :: group do p=1,size(config_name_phase) group = trim('current/constituent')//'/'//trim(config_name_phase(p)) - call HDF5_closeGroup(results_addGroup(group)) + call results_closeGroup(results_addGroup(group)) group = trim(group)//'/plastic' - call HDF5_closeGroup(results_addGroup(group)) + call results_closeGroup(results_addGroup(group)) select case(phase_plasticity(p)) case(PLASTICITY_ISOTROPIC_ID) diff --git a/src/plastic_disloUCLA.f90 b/src/constitutive_plastic_disloUCLA.f90 similarity index 95% rename from src/plastic_disloUCLA.f90 rename to src/constitutive_plastic_disloUCLA.f90 index 29ef0fcc9..a2b8f604d 100644 --- a/src/plastic_disloUCLA.f90 +++ b/src/constitutive_plastic_disloUCLA.f90 @@ -5,21 +5,9 @@ !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief crystal plasticity model for bcc metals, especially Tungsten !-------------------------------------------------------------------------------------------------- -module plastic_disloUCLA - use prec - use debug - use math - use IO - use material - use config - use lattice - use discretization - use results +submodule(constitutive) plastic_disloUCLA - implicit none - private - - real(pReal), parameter, private :: & + real(pReal), parameter :: & kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin enum, bind(c) @@ -33,14 +21,14 @@ module plastic_disloUCLA tau_pass_ID end enum - type, private :: tParameters + type :: tParameters real(pReal) :: & aTol_rho, & D, & !< grain size mu, & D_0, & !< prefactor for self-diffusion coefficient Q_cl !< activation energy for dislocation climb - real(pReal), dimension(:), allocatable :: & + real(pReal), dimension(:), allocatable :: & rho_mob_0, & !< initial dislocation density rho_dip_0, & !< initial dipole density b_sl, & !< magnitude of burgers vector [m] @@ -58,31 +46,31 @@ module plastic_disloUCLA kink_height, & !< height of the kink pair w, & !< width of the kink pair omega !< attempt frequency for kink pair nucleation - real(pReal), dimension(:,:), allocatable :: & + real(pReal), dimension(:,:), allocatable :: & h_sl_sl, & !< slip resistance from slip activity forestProjectionEdge - real(pReal), dimension(:,:,:), allocatable :: & + real(pReal), dimension(:,:,:), allocatable :: & Schmid, & nonSchmid_pos, & nonSchmid_neg integer :: & sum_N_sl !< total number of active slip system - integer, dimension(:), allocatable :: & + integer, dimension(:), allocatable :: & N_sl !< number of active slip systems for each family - integer(kind(undefined_ID)), dimension(:),allocatable :: & + integer(kind(undefined_ID)), dimension(:), allocatable :: & outputID !< ID of each post result output logical :: & dipoleFormation !< flag indicating consideration of dipole formation end type !< container type for internal constitutive parameters - type, private :: tDisloUCLAState + type :: tDisloUCLAState real(pReal), dimension(:,:), pointer :: & rho_mob, & rho_dip, & gamma_sl end type tDisloUCLAState - type, private :: tDisloUCLAdependentState + type :: tDisloUCLAdependentState real(pReal), dimension(:,:), allocatable :: & Lambda_sl, & threshold_stress @@ -90,20 +78,11 @@ module plastic_disloUCLA !-------------------------------------------------------------------------------------------------- ! containers for parameters and state - type(tParameters), allocatable, dimension(:), private :: param - type(tDisloUCLAState), allocatable, dimension(:), private :: & + type(tParameters), allocatable, dimension(:) :: param + type(tDisloUCLAState), allocatable, dimension(:) :: & dotState, & state - type(tDisloUCLAdependentState), allocatable, dimension(:), private :: dependentState - - public :: & - plastic_disloUCLA_init, & - plastic_disloUCLA_dependentState, & - plastic_disloUCLA_LpAndItsTangent, & - plastic_disloUCLA_dotState, & - plastic_disloUCLA_results - private :: & - kinetics + type(tDisloUCLAdependentState), allocatable, dimension(:) :: dependentState contains @@ -112,7 +91,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine plastic_disloUCLA_init() +module subroutine plastic_disloUCLA_init integer :: & Ninstance, & @@ -333,7 +312,7 @@ end subroutine plastic_disloUCLA_init !-------------------------------------------------------------------------------------------------- !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- -pure subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp, & +pure module subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp, & Mp,T,instance,of) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient @@ -376,7 +355,7 @@ end subroutine plastic_disloUCLA_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -subroutine plastic_disloUCLA_dotState(Mp,T,instance,of) +module subroutine plastic_disloUCLA_dotState(Mp,T,instance,of) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -436,7 +415,7 @@ end subroutine plastic_disloUCLA_dotState !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -subroutine plastic_disloUCLA_dependentState(instance,of) +module subroutine plastic_disloUCLA_dependentState(instance,of) integer, intent(in) :: & instance, & @@ -462,7 +441,7 @@ end subroutine plastic_disloUCLA_dependentState !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- -subroutine plastic_disloUCLA_results(instance,group) +module subroutine plastic_disloUCLA_results(instance,group) #if defined(PETSc) || defined(DAMASK_HDF5) integer, intent(in) :: instance @@ -615,4 +594,4 @@ pure subroutine kinetics(Mp,T,instance,of, & end subroutine kinetics -end module plastic_disloUCLA +end submodule plastic_disloUCLA diff --git a/src/plastic_dislotwin.f90 b/src/constitutive_plastic_dislotwin.f90 similarity index 98% rename from src/plastic_dislotwin.f90 rename to src/constitutive_plastic_dislotwin.f90 index 24b6bd3da..9f5dbb928 100644 --- a/src/plastic_dislotwin.f90 +++ b/src/constitutive_plastic_dislotwin.f90 @@ -7,19 +7,7 @@ !> @brief material subroutine incoprorating dislocation and twinning physics !> @details to be done !-------------------------------------------------------------------------------------------------- -module plastic_dislotwin - use prec - use debug - use math - use IO - use material - use config - use lattice - use discretization - use results - - implicit none - private +submodule(constitutive) plastic_dislotwin real(pReal), parameter :: & kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin @@ -156,14 +144,6 @@ module plastic_dislotwin state type(tDislotwinMicrostructure), allocatable, dimension(:) :: dependentState - public :: & - plastic_dislotwin_init, & - plastic_dislotwin_homogenizedC, & - plastic_dislotwin_dependentState, & - plastic_dislotwin_LpAndItsTangent, & - plastic_dislotwin_dotState, & - plastic_dislotwin_results - contains @@ -171,7 +151,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine plastic_dislotwin_init +module subroutine plastic_dislotwin_init integer :: & Ninstance, & @@ -576,14 +556,14 @@ end subroutine plastic_dislotwin_init !-------------------------------------------------------------------------------------------------- !> @brief returns the homogenized elasticity matrix !-------------------------------------------------------------------------------------------------- -function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC) +module function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC) real(pReal), dimension(6,6) :: & homogenizedC integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element integer :: i, & of @@ -615,7 +595,7 @@ end function plastic_dislotwin_homogenizedC !-------------------------------------------------------------------------------------------------- !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- -subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) +module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) real(pReal), dimension(3,3), intent(out) :: Lp real(pReal), dimension(3,3,3,3), intent(out) :: dLp_dMp @@ -730,7 +710,7 @@ end subroutine plastic_dislotwin_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -subroutine plastic_dislotwin_dotState(Mp,T,instance,of) +module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) real(pReal), dimension(3,3), intent(in):: & Mp !< Mandel stress @@ -833,7 +813,7 @@ end subroutine plastic_dislotwin_dotState !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -subroutine plastic_dislotwin_dependentState(T,instance,of) +module subroutine plastic_dislotwin_dependentState(T,instance,of) integer, intent(in) :: & instance, & @@ -925,11 +905,11 @@ end subroutine plastic_dislotwin_dependentState !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- -subroutine plastic_dislotwin_results(instance,group) +module subroutine plastic_dislotwin_results(instance,group) #if defined(PETSc) || defined(DAMASK_HDF5) - integer, intent(in) :: instance - character(len=*) :: group + integer, intent(in) :: instance + character(len=*), intent(in) :: group integer :: o associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) @@ -1183,4 +1163,4 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,of,& end subroutine kinetics_trans -end module plastic_dislotwin +end submodule plastic_dislotwin From 9882c3532a93221130f93ca6c8eb939de0dfe5bf Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 2 Dec 2019 22:26:08 +0100 Subject: [PATCH 04/20] avoid use of low-level HDF5 routines --- src/homogenization.f90 | 7 +++---- src/mesh_marc.f90 | 3 +-- 2 files changed, 4 insertions(+), 6 deletions(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 618c2faf8..0a61cbd3a 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -23,7 +23,6 @@ module homogenization use damage_local use damage_nonlocal use results - use HDF5_utilities implicit none private @@ -801,18 +800,18 @@ subroutine homogenization_results do p=1,size(config_name_homogenization) group = trim('current/materialpoint')//'/'//trim(config_name_homogenization(p)) - call HDF5_closeGroup(results_addGroup(group)) + call results_closeGroup(results_addGroup(group)) group = trim(group)//'/mech' - call HDF5_closeGroup(results_addGroup(group)) + call results_closeGroup(results_addGroup(group)) select case(material_homogenization_type(p)) case(HOMOGENIZATION_rgc_ID) call mech_RGC_results(homogenization_typeInstance(p),group) end select group = trim('current/materialpoint')//'/'//trim(config_name_homogenization(p))//'/generic' - call HDF5_closeGroup(results_addGroup(group)) + call results_closeGroup(results_addGroup(group)) !temp = reshape(materialpoint_F,[3,3,discretization_nIP*discretization_nElem]) !call results_writeDataset(group,temp,'F',& diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index f640baa72..b550f4229 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -18,7 +18,6 @@ module mesh use element use discretization use geometry_plastic_nonlocal - use HDF5_utilities use results implicit none @@ -152,7 +151,7 @@ subroutine writeGeometry(elemType, & #if defined(DAMASK_HDF5) call results_openJobFile - call HDF5_closeGroup(results_addGroup('geometry')) + call results_closeGroup(results_addGroup('geometry')) connectivity_temp = connectivity_elem call results_writeDataset('geometry',connectivity_temp,'T_e',& From ab1f0dc16b5d48576af8ee3721457f6ce84b99d3 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 4 Dec 2019 20:50:46 +0100 Subject: [PATCH 05/20] submodules allow inter-module communication --- src/constitutive.f90 | 71 ++++++++++++++++++- ....f90 => constitutive_plastic_nonlocal.f90} | 55 +++++--------- src/crystallite.f90 | 1 - src/lattice.f90 | 3 +- 4 files changed, 87 insertions(+), 43 deletions(-) rename src/{plastic_nonlocal.f90 => constitutive_plastic_nonlocal.f90} (98%) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 7f87d3cda..6ef5ee73d 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -6,6 +6,7 @@ module constitutive use prec use math + use rotations use debug use numerics use IO @@ -14,7 +15,6 @@ module constitutive use results use lattice use discretization - use plastic_nonlocal use geometry_plastic_nonlocal use source_thermal_dissipation use source_thermal_externalheat @@ -54,6 +54,9 @@ module constitutive module subroutine plastic_disloUCLA_init end subroutine plastic_disloUCLA_init + module subroutine plastic_nonlocal_init + end subroutine plastic_nonlocal_init + module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) real(pReal), dimension(3,3), intent(out) :: & @@ -123,6 +126,23 @@ module constitutive instance, & of end subroutine plastic_disloUCLA_LpAndItsTangent + + module subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dMp, & + Mp, Temperature, volume, 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, & + volume + integer, intent(in) :: & + ip, & !< current integration point + el !< current element number + end subroutine plastic_nonlocal_LpAndItsTangent module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of) @@ -182,6 +202,21 @@ module constitutive instance, & of end subroutine plastic_disloUCLA_dotState + + module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & + timestep,ip,el) + integer, intent(in) :: & + ip, & !< current integration point + el !< current element number + real(pReal), intent(in) :: & + Temperature, & !< temperature + timestep !< substepped crystallite time increment + real(pReal), dimension(3,3), intent(in) ::& + Mp !< MandelStress + real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: & + Fe, & !< elastic deformation gradient + Fp !< plastic deformation gradient + end subroutine plastic_nonlocal_dotState module subroutine plastic_dislotwin_dependentState(T,instance,of) @@ -197,6 +232,15 @@ module constitutive instance, & of end subroutine plastic_disloUCLA_dependentState + + module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) + integer, intent(in) :: & + ip, & + el + real(pReal), dimension(3,3), intent(in) :: & + Fe, & + Fp + end subroutine plastic_nonlocal_dependentState module subroutine plastic_kinehardening_deltaState(Mp,instance,of) @@ -206,6 +250,14 @@ module constitutive instance, & of end subroutine plastic_kinehardening_deltaState + + module subroutine plastic_nonlocal_deltaState(Mp,ip,el) + integer, intent(in) :: & + ip, & + el + real(pReal), dimension(3,3), intent(in) :: & + Mp + end subroutine plastic_nonlocal_deltaState module function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC) @@ -217,6 +269,14 @@ module constitutive el !< element end function plastic_dislotwin_homogenizedC + module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) + integer, intent(in) :: & + i, & + e + type(rotation), dimension(1,discretization_nIP,discretization_nElem), intent(in) :: & + orientation !< crystal orientation + end subroutine plastic_nonlocal_updateCompatibility + module subroutine plastic_isotropic_results(instance,group) integer, intent(in) :: instance @@ -242,10 +302,16 @@ module constitutive integer, intent(in) :: instance character(len=*), intent(in) :: group end subroutine plastic_disloUCLA_results + + module subroutine plastic_nonlocal_results(instance,group) + integer, intent(in) :: instance + character(len=*), intent(in) :: group + end subroutine plastic_nonlocal_results end interface public :: & + plastic_nonlocal_updateCompatibility, & constitutive_init, & constitutive_homogenizedC, & constitutive_microstructure, & @@ -446,8 +512,7 @@ end subroutine constitutive_microstructure ! Mp in, dLp_dMp out !-------------------------------------------------------------------------------------------------- subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & - S, Fi, ipc, ip, el) - + S, Fi, ipc, ip, el) integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point diff --git a/src/plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 similarity index 98% rename from src/plastic_nonlocal.f90 rename to src/constitutive_plastic_nonlocal.f90 index 5a3f8c173..099a1304e 100644 --- a/src/plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -4,17 +4,8 @@ !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @brief material subroutine for plasticity including dislocation flux !-------------------------------------------------------------------------------------------------- -module plastic_nonlocal - use prec - use IO - use math - use debug - use material - use lattice +submodule(constitutive) plastic_nonlocal use rotations - use config - use lattice - use discretization use geometry_plastic_nonlocal, only: & nIPneighbors => geometry_plastic_nonlocal_nIPneighbors, & IPneighborhood => geometry_plastic_nonlocal_IPneighborhood, & @@ -22,8 +13,7 @@ module plastic_nonlocal IParea => geometry_plastic_nonlocal_IParea0, & IPareaNormal => geometry_plastic_nonlocal_IPareaNormal0 - implicit none - private + real(pReal), parameter :: & KB = 1.38e-23_pReal !< Physical parameter, Boltzmann constant in J/Kelvin @@ -62,7 +52,7 @@ module plastic_nonlocal !END DEPRECATED real(pReal), dimension(:,:,:,:,:,:), allocatable :: & - compatibility !< slip system compatibility between me and my neighbors + compatibility !< slip system compatibility between me and my neighbors enum, bind(c) enumerator :: & @@ -148,7 +138,7 @@ module plastic_nonlocal nonSchmid_neg !< combined projection of Schmid and non-Schmid contributions to the resolved shear stress (only for screws) integer :: & totalNslip - integer, dimension(:) ,allocatable :: & + integer, dimension(:) ,allocatable:: & Nslip,& colinearSystem !< colinear system to the active slip system (only valid for fcc!) @@ -204,18 +194,7 @@ module plastic_nonlocal integer(kind(undefined_ID)), dimension(:,:), allocatable :: & plastic_nonlocal_outputID !< ID of each post result output - - public :: & - plastic_nonlocal_init, & - plastic_nonlocal_dependentState, & - plastic_nonlocal_LpAndItsTangent, & - plastic_nonlocal_dotState, & - plastic_nonlocal_deltaState, & - plastic_nonlocal_updateCompatibility, & - plastic_nonlocal_results - - private :: & - plastic_nonlocal_kinetics + contains @@ -223,7 +202,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine plastic_nonlocal_init +module subroutine plastic_nonlocal_init character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] integer, dimension(0), parameter :: emptyIntArray = [integer::] @@ -752,7 +731,7 @@ end subroutine plastic_nonlocal_init !-------------------------------------------------------------------------------------------------- !> @brief calculates quantities characterizing the microstructure !-------------------------------------------------------------------------------------------------- -subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) +module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) integer, intent(in) :: & ip, & @@ -1114,7 +1093,7 @@ end subroutine plastic_nonlocal_kinetics !-------------------------------------------------------------------------------------------------- !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- -subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dMp, & +module subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dMp, & Mp, Temperature, volume, ip, el) integer, intent(in) :: & @@ -1244,7 +1223,7 @@ end subroutine plastic_nonlocal_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief (instantaneous) incremental change of microstructure !-------------------------------------------------------------------------------------------------- -subroutine plastic_nonlocal_deltaState(Mp,ip,el) +module subroutine plastic_nonlocal_deltaState(Mp,ip,el) integer, intent(in) :: & ip, & @@ -1360,8 +1339,8 @@ end subroutine plastic_nonlocal_deltaState !--------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !--------------------------------------------------------------------------------------------------- -subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & - timestep,ip,el) +module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & + timestep,ip,el) integer, intent(in) :: & ip, & !< current integration point @@ -1809,13 +1788,13 @@ end subroutine plastic_nonlocal_dotState ! plane normals and signed cosine of the angle between the slip directions. Only the largest values ! that sum up to a total of 1 are considered, all others are set to zero. !-------------------------------------------------------------------------------------------------- -subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) +module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) integer, intent(in) :: & i, & e type(rotation), dimension(1,discretization_nIP,discretization_nElem), intent(in) :: & - orientation ! crystal orientation in quaternions + orientation ! crystal orientation integer :: & Nneighbors, & ! number of neighbors @@ -1974,13 +1953,13 @@ end function getRho !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- -subroutine plastic_nonlocal_results(instance,group) +module subroutine plastic_nonlocal_results(instance,group) #if defined(PETSc) || defined(DAMASK_HDF5) use results, only: & results_writeDataset - integer, intent(in) :: instance - character(len=*) :: group + integer, intent(in) :: instance + character(len=*),intent(in) :: group integer :: o associate(prm => param(instance),dst => microstructure(instance),stt=>state(instance)) @@ -2047,4 +2026,4 @@ subroutine plastic_nonlocal_results(instance,group) end subroutine plastic_nonlocal_results -end module plastic_nonlocal +end submodule plastic_nonlocal diff --git a/src/crystallite.f90 b/src/crystallite.f90 index a17213f9e..c18588218 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -21,7 +21,6 @@ module crystallite use constitutive use discretization use lattice - use plastic_nonlocal use results implicit none diff --git a/src/lattice.f90 b/src/lattice.f90 index fada61392..dfeac6ce1 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -483,7 +483,8 @@ module lattice lattice_slip_normal, & lattice_slip_direction, & lattice_slip_transverse, & - lattice_labels_slip + lattice_labels_slip, & + lattice_labels_twin contains !-------------------------------------------------------------------------------------------------- From e70c56701fbdbc9a0f57461a66ea036500f9ccb8 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 13 Jan 2020 21:09:25 +0100 Subject: [PATCH 06/20] not needed --- src/constitutive_plastic_nonlocal.f90 | 1 - 1 file changed, 1 deletion(-) diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index ce052b050..42b8a34f5 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -5,7 +5,6 @@ !> @brief material subroutine for plasticity including dislocation flux !-------------------------------------------------------------------------------------------------- submodule(constitutive) plastic_nonlocal - use rotations use geometry_plastic_nonlocal, only: & nIPneighbors => geometry_plastic_nonlocal_nIPneighbors, & IPneighborhood => geometry_plastic_nonlocal_IPneighborhood, & From 27bc23c2e11db139b028728c92dcd7f4ae27bf46 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 14 Jan 2020 07:55:18 +0100 Subject: [PATCH 07/20] missing renames wondering how this has passed the syntax check earlier --- src/commercialFEM_fileList.f90 | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index 31bddd141..c22a37a0f 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -36,14 +36,14 @@ #include "kinematics_cleavage_opening.f90" #include "kinematics_slipplane_opening.f90" #include "kinematics_thermal_expansion.f90" -#include "plastic_none.f90" -#include "plastic_isotropic.f90" -#include "plastic_phenopowerlaw.f90" -#include "plastic_kinematichardening.f90" -#include "plastic_dislotwin.f90" -#include "plastic_disloUCLA.f90" -#include "plastic_nonlocal.f90" #include "constitutive.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_disloUCLA.f90" +#include "constitutive_plastic_nonlocal.f90" #include "crystallite.f90" #include "thermal_isothermal.f90" #include "thermal_adiabatic.f90" From c8835947dbddab826945e76187a617245eb6b00c Mon Sep 17 00:00:00 2001 From: Test User Date: Fri, 24 Jan 2020 22:16:33 +0100 Subject: [PATCH 08/20] [skip ci] updated version information after successful test of v2.0.3-1545-ge63c9b94 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 3339b66fc..8dda7a7de 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-1531-g3d6ec695 +v2.0.3-1545-ge63c9b94 From 2100742a31a93365adb2c4adf4b4b7d0d3602179 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 26 Jan 2020 08:17:25 +0100 Subject: [PATCH 09/20] better have no support then untested support --- src/mesh_marc.f90 | 27 ++++++++++++--------------- 1 file changed, 12 insertions(+), 15 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 22f978609..869677b19 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -561,26 +561,23 @@ subroutine inputRead_elemType(elem, & select case (IO_lc(what)) case ( '6') mapElemtype = 1 ! Two-dimensional Plane Strain Triangle - case ( '155', & - '125', & - '128') + case ( '125') ! 155, 128 (need test) mapElemtype = 2 ! Two-dimensional Plane Strain triangle (155: cubic shape function, 125/128: second order isoparametric) - case ( '11') - mapElemtype = 3 ! Arbitrary Quadrilateral Plane-strain + !case ( '11') ! need test + ! mapElemtype = 3 ! Arbitrary Quadrilateral Plane-strain case ( '27') mapElemtype = 4 ! Plane Strain, Eight-node Distorted Quadrilateral case ( '54') mapElemtype = 5 ! Plane Strain, Eight-node Distorted Quadrilateral with reduced integration - case ( '134') - mapElemtype = 6 ! Three-dimensional Four-node Tetrahedron - case ( '157') - mapElemtype = 7 ! Three-dimensional, Low-order, Tetrahedron, Herrmann Formulations - case ( '127') - mapElemtype = 8 ! Three-dimensional Ten-node Tetrahedron - case ( '136') - mapElemtype = 9 ! Three-dimensional Arbitrarily Distorted Pentahedral - case ( '117', & - '123') + !case ( '134') ! need test + ! mapElemtype = 6 ! Three-dimensional Four-node Tetrahedron + !case ( '157') ! need test + ! mapElemtype = 7 ! Three-dimensional, Low-order, Tetrahedron, Herrmann Formulations + !case ( '127') ! need test + ! mapElemtype = 8 ! Three-dimensional Ten-node Tetrahedron + !case ( '136') ! need test + ! mapElemtype = 9 ! Three-dimensional Arbitrarily Distorted Pentahedral + case ( '117') ! 123 (need test) mapElemtype = 10 ! Three-dimensional Arbitrarily Distorted linear hexahedral with reduced integration case ( '7') mapElemtype = 11 ! Three-dimensional Arbitrarily Distorted Brick From 7671e257bd489704a015fb7d1f9b568691fdfa04 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 26 Jan 2020 08:17:59 +0100 Subject: [PATCH 10/20] reorder elem 11 (Marc 7) to match 12 (Marc 57) --- src/element.f90 | 92 +++++++++++++++++++++++++------------------------ 1 file changed, 47 insertions(+), 45 deletions(-) diff --git a/src/element.f90 b/src/element.f90 index a38036f93..023259c9b 100644 --- a/src/element.f90 +++ b/src/element.f90 @@ -482,33 +482,34 @@ module element integer, dimension(NNODE(11),NCELLNODE(GEOMTYPE(11))), parameter :: CELLNODEPARENTNODEWEIGHTS11 = & reshape([& - 1, 0, 0, 0, 0, 0, 0, 0, & ! - 0, 1, 0, 0, 0, 0, 0, 0, & ! - 0, 0, 1, 0, 0, 0, 0, 0, & ! - 0, 0, 0, 1, 0, 0, 0, 0, & ! + 1, 0, 0, 0, 0, 0, 0, 0, & ! 1 + 0, 1, 0, 0, 0, 0, 0, 0, & ! 2 + 0, 0, 1, 0, 0, 0, 0, 0, & ! 3 + 0, 0, 0, 1, 0, 0, 0, 0, & ! 4 0, 0, 0, 0, 1, 0, 0, 0, & ! 5 - 0, 0, 0, 0, 0, 1, 0, 0, & ! - 0, 0, 0, 0, 0, 0, 1, 0, & ! - 0, 0, 0, 0, 0, 0, 0, 1, & ! - 1, 1, 0, 0, 0, 0, 0, 0, & ! + 0, 0, 0, 0, 0, 1, 0, 0, & ! 6 + 0, 0, 0, 0, 0, 0, 1, 0, & ! 7 + 0, 0, 0, 0, 0, 0, 0, 1, & ! 8 + 1, 1, 0, 0, 0, 0, 0, 0, & ! 9 0, 1, 1, 0, 0, 0, 0, 0, & ! 10 - 0, 0, 1, 1, 0, 0, 0, 0, & ! - 1, 0, 0, 1, 0, 0, 0, 0, & ! - 1, 0, 0, 0, 1, 0, 0, 0, & ! - 0, 1, 0, 0, 0, 1, 0, 0, & ! - 0, 0, 1, 0, 0, 0, 1, 0, & ! 15 - 0, 0, 0, 1, 0, 0, 0, 1, & ! - 0, 0, 0, 0, 1, 1, 0, 0, & ! - 0, 0, 0, 0, 0, 1, 1, 0, & ! - 0, 0, 0, 0, 0, 0, 1, 1, & ! + 0, 0, 1, 1, 0, 0, 0, 0, & ! 11 + 1, 0, 0, 1, 0, 0, 0, 0, & ! 12 + 0, 0, 0, 0, 1, 1, 0, 0, & ! 17 + 0, 0, 0, 0, 0, 1, 1, 0, & ! 18 + 0, 0, 0, 0, 0, 0, 1, 1, & ! 19 0, 0, 0, 0, 1, 0, 0, 1, & ! 20 - 1, 1, 1, 1, 0, 0, 0, 0, & ! - 1, 1, 0, 0, 1, 1, 0, 0, & ! - 0, 1, 1, 0, 0, 1, 1, 0, & ! - 0, 0, 1, 1, 0, 0, 1, 1, & ! + 1, 0, 0, 0, 1, 0, 0, 0, & ! 13 + 0, 1, 0, 0, 0, 1, 0, 0, & ! 14 + 0, 0, 1, 0, 0, 0, 1, 0, & ! 15 + 0, 0, 0, 1, 0, 0, 0, 1, & ! 16 +! + 1, 1, 1, 1, 0, 0, 0, 0, & ! 21 + 1, 1, 0, 0, 1, 1, 0, 0, & ! 22 + 0, 1, 1, 0, 0, 1, 1, 0, & ! 23 + 0, 0, 1, 1, 0, 0, 1, 1, & ! 24 1, 0, 0, 1, 1, 0, 0, 1, & ! 25 - 0, 0, 0, 0, 1, 1, 1, 1, & ! - 1, 1, 1, 1, 1, 1, 1, 1 & ! + 0, 0, 0, 0, 1, 1, 1, 1, & ! 26 + 1, 1, 1, 1, 1, 1, 1, 1 & ! 27 #if !defined(__GFORTRAN__) ],shape(CELLNODEPARENTNODEWEIGHTS11)) !< 3D 8node 8ip #else @@ -517,33 +518,34 @@ module element integer, dimension(NNODE(12),NCELLNODE(GEOMTYPE(12))), parameter :: CELLNODEPARENTNODEWEIGHTS12 = & reshape([& - 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! + 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! 1 + 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! 2 + 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! 3 + 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! 4 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! 5 - 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! + 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! 6 + 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! 7 + 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! 8 + 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! 9 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! 10 - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, & ! - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, & ! + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! 11 + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, & ! 12 + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, & ! 13 + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, & ! 14 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, & ! 15 - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, & ! - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, & ! - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, & ! - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, & ! + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, & ! 16 + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, & ! 17 + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, & ! 18 + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, & ! 19 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, & ! 20 - 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 1, 1, 0, 0, 1, 1, 0, 0, 2, 0, 0, 0, 2, 0, 0, 0, 2, 2, 0, 0, & ! - 0, 1, 1, 0, 0, 1, 1, 0, 0, 2, 0, 0, 0, 2, 0, 0, 0, 2, 2, 0, & ! - 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 2, 0, 0, 0, 2, 0, 0, 0, 2, 2, & ! +! + 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, & ! 21 + 1, 1, 0, 0, 1, 1, 0, 0, 2, 0, 0, 0, 2, 0, 0, 0, 2, 2, 0, 0, & ! 22 + 0, 1, 1, 0, 0, 1, 1, 0, 0, 2, 0, 0, 0, 2, 0, 0, 0, 2, 2, 0, & ! 23 + 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 2, 0, 0, 0, 2, 0, 0, 0, 2, 2, & ! 24 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 2, 0, 0, 0, 2, 2, 0, 0, 2, & ! 25 - 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 2, 2, 0, 0, 0, 0, & ! - 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 & ! + 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 2, 2, 0, 0, 0, 0, & ! 26 + 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 & ! 27 #if !defined(__GFORTRAN__) ],shape(CELLNODEPARENTNODEWEIGHTS12)) !< 3D 20node 8ip #else From b16f4155d487741c20088bb5256f7cf50e05650d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 26 Jan 2020 09:06:45 +0100 Subject: [PATCH 11/20] cell definition follows correct order for 11/12 (7/57 in MSC.Marc) --- src/element.f90 | 64 ++++++++++++++++++++++++------------------------- 1 file changed, 32 insertions(+), 32 deletions(-) diff --git a/src/element.f90 b/src/element.f90 index 023259c9b..91b687c2e 100644 --- a/src/element.f90 +++ b/src/element.f90 @@ -482,34 +482,34 @@ module element integer, dimension(NNODE(11),NCELLNODE(GEOMTYPE(11))), parameter :: CELLNODEPARENTNODEWEIGHTS11 = & reshape([& - 1, 0, 0, 0, 0, 0, 0, 0, & ! 1 - 0, 1, 0, 0, 0, 0, 0, 0, & ! 2 - 0, 0, 1, 0, 0, 0, 0, 0, & ! 3 - 0, 0, 0, 1, 0, 0, 0, 0, & ! 4 + 1, 0, 0, 0, 0, 0, 0, 0, & ! + 0, 1, 0, 0, 0, 0, 0, 0, & ! + 0, 0, 1, 0, 0, 0, 0, 0, & ! + 0, 0, 0, 1, 0, 0, 0, 0, & ! 0, 0, 0, 0, 1, 0, 0, 0, & ! 5 - 0, 0, 0, 0, 0, 1, 0, 0, & ! 6 - 0, 0, 0, 0, 0, 0, 1, 0, & ! 7 - 0, 0, 0, 0, 0, 0, 0, 1, & ! 8 - 1, 1, 0, 0, 0, 0, 0, 0, & ! 9 + 0, 0, 0, 0, 0, 1, 0, 0, & ! + 0, 0, 0, 0, 0, 0, 1, 0, & ! + 0, 0, 0, 0, 0, 0, 0, 1, & ! + 1, 1, 0, 0, 0, 0, 0, 0, & ! 0, 1, 1, 0, 0, 0, 0, 0, & ! 10 - 0, 0, 1, 1, 0, 0, 0, 0, & ! 11 - 1, 0, 0, 1, 0, 0, 0, 0, & ! 12 - 0, 0, 0, 0, 1, 1, 0, 0, & ! 17 - 0, 0, 0, 0, 0, 1, 1, 0, & ! 18 - 0, 0, 0, 0, 0, 0, 1, 1, & ! 19 - 0, 0, 0, 0, 1, 0, 0, 1, & ! 20 - 1, 0, 0, 0, 1, 0, 0, 0, & ! 13 - 0, 1, 0, 0, 0, 1, 0, 0, & ! 14 - 0, 0, 1, 0, 0, 0, 1, 0, & ! 15 - 0, 0, 0, 1, 0, 0, 0, 1, & ! 16 + 0, 0, 1, 1, 0, 0, 0, 0, & ! + 1, 0, 0, 1, 0, 0, 0, 0, & ! + 0, 0, 0, 0, 1, 1, 0, 0, & ! 17 => 13 + 0, 0, 0, 0, 0, 1, 1, 0, & ! 18 => 14 + 0, 0, 0, 0, 0, 0, 1, 1, & ! 19 => 15 + 0, 0, 0, 0, 1, 0, 0, 1, & ! 20 => 16 + 1, 0, 0, 0, 1, 0, 0, 0, & ! 13 => 17 + 0, 1, 0, 0, 0, 1, 0, 0, & ! 14 => 18 + 0, 0, 1, 0, 0, 0, 1, 0, & ! 15 => 19 + 0, 0, 0, 1, 0, 0, 0, 1, & ! 16 => 20 ! - 1, 1, 1, 1, 0, 0, 0, 0, & ! 21 - 1, 1, 0, 0, 1, 1, 0, 0, & ! 22 - 0, 1, 1, 0, 0, 1, 1, 0, & ! 23 - 0, 0, 1, 1, 0, 0, 1, 1, & ! 24 + 1, 1, 1, 1, 0, 0, 0, 0, & ! + 1, 1, 0, 0, 1, 1, 0, 0, & ! + 0, 1, 1, 0, 0, 1, 1, 0, & ! + 0, 0, 1, 1, 0, 0, 1, 1, & ! 1, 0, 0, 1, 1, 0, 0, 1, & ! 25 - 0, 0, 0, 0, 1, 1, 1, 1, & ! 26 - 1, 1, 1, 1, 1, 1, 1, 1 & ! 27 + 0, 0, 0, 0, 1, 1, 1, 1, & ! + 1, 1, 1, 1, 1, 1, 1, 1 & ! #if !defined(__GFORTRAN__) ],shape(CELLNODEPARENTNODEWEIGHTS11)) !< 3D 8node 8ip #else @@ -720,14 +720,14 @@ module element integer, dimension(NCELLNODEPERCELL(CELLTYPE(9)),NIP(9)), parameter :: CELL9 = & reshape([& - 1, 9,21,12,13,22,27,25, & - 9, 2,10,21,22,14,23,27, & - 12,21,11, 4,25,27,24,16, & - 21,10, 3,11,27,23,15,24, & - 13,22,27,25, 5,17,26,20, & - 22,14,23,27,17, 6,18,26, & - 25,27,24,16,20,26,19, 8, & - 27,23,15,24,26,18, 7,19 & + 1, 9,21,12,17,22,27,25, & + 9, 2,10,21,22,18,23,27, & + 12,21,11, 4,25,27,24,20, & + 21,10, 3,11,27,23,19,24, & + 17,22,27,25, 5,13,26,16, & + 22,18,23,27,13, 6,14,26, & + 25,27,24,20,16,26,15, 8, & + 27,23,19,24,26,14, 7,15 & #if !defined(__GFORTRAN__) ],shape(CELL9)) #else From ffea69955e7a674e4be3b38b4249086debfa6b67 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 26 Jan 2020 09:14:51 +0100 Subject: [PATCH 12/20] polishing and testing element 11/12 (7/57 for MSC.Marc) have the same geometry type, i.e. both are a hexahedral with 8 integration points/cells Type 11 is linear (8 nodes), type 12 quadradic with reduced integration (20 nodes). The modified definition ensures that the cell nodes 9-20 of element 11 are in the same order as the real nodes 9-20 of element 12. Real nodes 1-8 (corners) and 21-27 needed no modification. Notes: * Documentation on https://damask.mpie.de/Documentation/ElementType is now outdated. * Element defition in MSC.Marc manual volume B (2001 version) is confusing because element numbering is sometimes clowise and sometimes counterclockwise. The latter one seems to be correct --- PRIVATE | 2 +- src/element.f90 | 62 ++++++++++++++++++++++++------------------------- 2 files changed, 31 insertions(+), 33 deletions(-) diff --git a/PRIVATE b/PRIVATE index 5be24da73..11aea4aa2 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 5be24da737ceb89599c7b995c605cab08fdc05c2 +Subproject commit 11aea4aa272ef037c189a158c4f7b2af923fdf30 diff --git a/src/element.f90 b/src/element.f90 index 91b687c2e..1159d401a 100644 --- a/src/element.f90 +++ b/src/element.f90 @@ -494,15 +494,14 @@ module element 0, 1, 1, 0, 0, 0, 0, 0, & ! 10 0, 0, 1, 1, 0, 0, 0, 0, & ! 1, 0, 0, 1, 0, 0, 0, 0, & ! - 0, 0, 0, 0, 1, 1, 0, 0, & ! 17 => 13 - 0, 0, 0, 0, 0, 1, 1, 0, & ! 18 => 14 - 0, 0, 0, 0, 0, 0, 1, 1, & ! 19 => 15 - 0, 0, 0, 0, 1, 0, 0, 1, & ! 20 => 16 - 1, 0, 0, 0, 1, 0, 0, 0, & ! 13 => 17 - 0, 1, 0, 0, 0, 1, 0, 0, & ! 14 => 18 - 0, 0, 1, 0, 0, 0, 1, 0, & ! 15 => 19 - 0, 0, 0, 1, 0, 0, 0, 1, & ! 16 => 20 -! + 0, 0, 0, 0, 1, 1, 0, 0, & ! + 0, 0, 0, 0, 0, 1, 1, 0, & ! + 0, 0, 0, 0, 0, 0, 1, 1, & ! 15 + 0, 0, 0, 0, 1, 0, 0, 1, & ! + 1, 0, 0, 0, 1, 0, 0, 0, & ! + 0, 1, 0, 0, 0, 1, 0, 0, & ! + 0, 0, 1, 0, 0, 0, 1, 0, & ! + 0, 0, 0, 1, 0, 0, 0, 1, & ! 20 1, 1, 1, 1, 0, 0, 0, 0, & ! 1, 1, 0, 0, 1, 1, 0, 0, & ! 0, 1, 1, 0, 0, 1, 1, 0, & ! @@ -518,34 +517,33 @@ module element integer, dimension(NNODE(12),NCELLNODE(GEOMTYPE(12))), parameter :: CELLNODEPARENTNODEWEIGHTS12 = & reshape([& - 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! 1 - 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! 2 - 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! 3 - 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! 4 + 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! + 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! + 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! + 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! 5 - 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! 6 - 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! 7 - 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! 8 - 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! 9 + 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! + 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! + 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! + 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! 10 - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! 11 - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, & ! 12 - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, & ! 13 - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, & ! 14 + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, & ! + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, & ! + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, & ! 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, & ! 15 - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, & ! 16 - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, & ! 17 - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, & ! 18 - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, & ! 19 + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, & ! + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, & ! + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, & ! + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, & ! 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, & ! 20 -! - 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, & ! 21 - 1, 1, 0, 0, 1, 1, 0, 0, 2, 0, 0, 0, 2, 0, 0, 0, 2, 2, 0, 0, & ! 22 - 0, 1, 1, 0, 0, 1, 1, 0, 0, 2, 0, 0, 0, 2, 0, 0, 0, 2, 2, 0, & ! 23 - 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 2, 0, 0, 0, 2, 0, 0, 0, 2, 2, & ! 24 + 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, & ! + 1, 1, 0, 0, 1, 1, 0, 0, 2, 0, 0, 0, 2, 0, 0, 0, 2, 2, 0, 0, & ! + 0, 1, 1, 0, 0, 1, 1, 0, 0, 2, 0, 0, 0, 2, 0, 0, 0, 2, 2, 0, & ! + 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 2, 0, 0, 0, 2, 0, 0, 0, 2, 2, & ! 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 2, 0, 0, 0, 2, 2, 0, 0, 2, & ! 25 - 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 2, 2, 0, 0, 0, 0, & ! 26 - 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 & ! 27 + 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 2, 2, 0, 0, 0, 0, & ! + 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 & ! #if !defined(__GFORTRAN__) ],shape(CELLNODEPARENTNODEWEIGHTS12)) !< 3D 20node 8ip #else From 0eba4e39ccfc6dd99f9bb939c38bbad226a88c06 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 26 Jan 2020 14:18:16 +0100 Subject: [PATCH 13/20] trying to make the whole procedure understandable --- src/element.f90 | 113 ++++++++++++++++++++++++------------------------ 1 file changed, 57 insertions(+), 56 deletions(-) diff --git a/src/element.f90 b/src/element.f90 index 1159d401a..6ce23a9da 100644 --- a/src/element.f90 +++ b/src/element.f90 @@ -43,38 +43,39 @@ module element integer, dimension(NELEMTYPE), parameter :: NNODE = & [ & - 3, & ! 2D 3node 1ip - 6, & ! 2D 6node 3ip - 4, & ! 2D 4node 4ip - 8, & ! 2D 8node 9ip - 8, & ! 2D 8node 4ip - !-------------------- - 4, & ! 3D 4node 1ip - 5, & ! 3D 5node 4ip - 10, & ! 3D 10node 4ip - 6, & ! 3D 6node 6ip - 8, & ! 3D 8node 1ip - 8, & ! 3D 8node 8ip - 20, & ! 3D 20node 8ip - 20 & ! 3D 20node 27ip + 3, & ! 2D, 1 IP + 6, & ! 2D, 3 IP + 4, & ! 2D, 4 IP + 8, & ! 2D, 9 IP + 8, & ! 2D, 4 IP + !---------------------- + 4, & ! 3D, 1 IP + 5, & ! 3D, 4 IP + 10, & ! 3D, 4 IP + 6, & ! 3D, 6 IP + 8, & ! 3D, 1 IP + 8, & ! 3D, 8 IP + 20, & ! 3D, 8 IP + 20 & ! 3D, 27 IP ] !< number of nodes that constitute a specific type of element integer, dimension(NELEMTYPE), parameter :: GEOMTYPE = & [ & - 1, & - 2, & - 3, & - 4, & - 3, & - 5, & - 6, & - 6, & - 7, & - 8, & - 9, & - 9, & - 10 & - ] !< geometry type of particular element type + 1, & ! 1 triangle + 2, & ! 3 quadrilaterals + 3, & ! 4 quadrilaterals + 4, & ! 9 quadrilaterals + 3, & ! 4 quadrilaterals + !---------------------- + 5, & ! 1 tetrahedron + 6, & ! 4 hexahedrons + 6, & ! 4 hexahedrons + 7, & ! 6 hexahedrons + 8, & ! 1 hexahedron + 9, & ! 8 hexahedrons + 9, & ! 8 hexahedrons + 10 & ! 27 hexahedrons + ] !< geometry type (same number of cell nodes and IPs) integer, dimension(maxval(GEOMTYPE)), parameter :: NCELLNODE = & [ & @@ -88,7 +89,7 @@ module element 8, & 27, & 64 & - ] !< number of cell nodes in a specific geometry type + ] !< number of cell nodes integer, dimension(maxval(GEOMTYPE)), parameter :: NIP = & [ & @@ -102,45 +103,45 @@ module element 1, & 8, & 27 & - ] !< number of IPs in a specific geometry type + ] !< number of IPs integer, dimension(maxval(GEOMTYPE)), parameter :: CELLTYPE = & [ & - 1, & ! 2D 3node - 2, & ! 2D 4node - 2, & ! 2D 4node - 2, & ! 2D 4node - 3, & ! 3D 4node - 4, & ! 3D 8node - 4, & ! 3D 8node - 4, & ! 3D 8node - 4, & ! 3D 8node - 4 & ! 3D 8node - ] !< cell type that is used by each geometry type + 1, & ! 2D, 3 node (Triangle) + 2, & ! 2D, 4 node (Quadrilateral) + 2, & ! - " - + 2, & ! - " - + 3, & ! 3D, 4 node (Tetrahedron) + 4, & ! 3D, 4 node (Hexahedron) + 4, & ! - " - + 4, & ! - " - + 4, & ! - " - + 4 & ! - " - + ] !< cell type integer, dimension(maxval(CELLTYPE)), parameter :: NIPNEIGHBOR = & [ & - 3, & ! 2D 3node - 4, & ! 2D 4node - 4, & ! 3D 4node - 6 & ! 3D 8node - ] !< number of ip neighbors / cell faces in a specific cell type + 3, & + 4, & + 4, & + 6 & + ] !< number of ip neighbors / cell faces integer, dimension(maxval(CELLTYPE)), parameter :: NCELLNODEPERCELLFACE = & [ & - 2, & ! 2D 3node - 2, & ! 2D 4node - 3, & ! 3D 4node - 4 & ! 3D 8node - ] !< number of cell nodes in a specific cell type + 2, & + 2, & + 3, & + 4 & + ] !< number of cell nodes per face integer, dimension(maxval(CELLTYPE)), parameter :: NCELLNODEPERCELL = & [ & - 3, & ! 2D 3node - 4, & ! 2D 4node - 4, & ! 3D 4node - 8 & ! 3D 8node - ] !< number of cell nodes in a specific cell type + 3, & + 4, & + 4, & + 8 & + ] !< number of total cell nodes ! *** IPneighbor *** ! list of the neighborhood of each IP. From 8344e23c65e40ea80d385cae185f7db63886ccab Mon Sep 17 00:00:00 2001 From: Test User Date: Mon, 27 Jan 2020 16:54:01 +0100 Subject: [PATCH 14/20] [skip ci] updated version information after successful test of v2.0.3-1601-gc433e244 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 8dda7a7de..66536aad6 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-1545-ge63c9b94 +v2.0.3-1601-gc433e244 From 55e53536f2510887ae9aded05ce74a723f434e66 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 31 Jan 2020 21:37:18 +0100 Subject: [PATCH 15/20] fixing indentation always 2 spaces, not 1 for the first level --- src/constitutive_plastic_dislotwin.f90 | 1933 ++++++++++---------- src/constitutive_plastic_kinehardening.f90 | 792 ++++---- src/constitutive_plastic_none.f90 | 36 +- src/constitutive_plastic_phenopowerlaw.f90 | 930 +++++----- 4 files changed, 1846 insertions(+), 1845 deletions(-) diff --git a/src/constitutive_plastic_dislotwin.f90 b/src/constitutive_plastic_dislotwin.f90 index 0972d9a2b..068323f3e 100644 --- a/src/constitutive_plastic_dislotwin.f90 +++ b/src/constitutive_plastic_dislotwin.f90 @@ -8,141 +8,141 @@ !> @details to be done !-------------------------------------------------------------------------------------------------- submodule(constitutive) plastic_dislotwin - - real(pReal), parameter :: & - kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin - - enum, bind(c) - enumerator :: & - undefined_ID, & - rho_mob_ID, & - rho_dip_ID, & - dot_gamma_sl_ID, & - gamma_sl_ID, & - Lambda_sl_ID, & - resolved_stress_slip_ID, & - tau_pass_ID, & - edge_dipole_distance_ID, & - f_tw_ID, & - Lambda_tw_ID, & - resolved_stress_twin_ID, & - tau_hat_tw_ID, & - f_tr_ID - end enum - - type :: tParameters - real(pReal) :: & - mu, & - nu, & - D0, & !< prefactor for self-diffusion coefficient - Qsd, & !< activation energy for dislocation climb - omega, & !< frequency factor for dislocation climb - D, & !< grain size - p_sb, & !< p-exponent in shear band velocity - q_sb, & !< q-exponent in shear band velocity - CEdgeDipMinDistance, & !< - i_tw, & !< - tau_0, & !< strength due to elements in solid solution - L_tw, & !< Length of twin nuclei in Burgers vectors - L_tr, & !< Length of trans nuclei in Burgers vectors - xc_twin, & !< critical distance for formation of twin nucleus - xc_trans, & !< critical distance for formation of trans nucleus - V_cs, & !< cross slip volume - sbResistance, & !< value for shearband resistance (might become an internal state variable at some point) - sbVelocity, & !< value for shearband velocity_0 - sbQedge, & !< activation energy for shear bands - SFE_0K, & !< stacking fault energy at zero K - dSFE_dT, & !< temperature dependance of stacking fault energy - aTol_rho, & !< absolute tolerance for integration of dislocation density - aTol_f_tw, & !< absolute tolerance for integration of twin volume fraction - aTol_f_tr, & !< absolute tolerance for integration of trans volume fraction - gamma_fcc_hex, & !< Free energy difference between austensite and martensite - i_tr, & !< - h !< Stack height of hex nucleus - real(pReal), dimension(:), allocatable :: & - rho_mob_0, & !< initial unipolar dislocation density per slip system - rho_dip_0, & !< initial dipole dislocation density per slip system - b_sl, & !< absolute length of burgers vector [m] for each slip system - b_tw, & !< absolute length of burgers vector [m] for each twin system - b_tr, & !< absolute length of burgers vector [m] for each transformation system - Delta_F,& !< activation energy for glide [J] for each slip system - v0, & !< dislocation velocity prefactor [m/s] for each slip system - dot_N_0_tw, & !< twin nucleation rate [1/m³s] for each twin system - dot_N_0_tr, & !< trans nucleation rate [1/m³s] for each trans system - t_tw, & !< twin thickness [m] for each twin system - CLambdaSlip, & !< Adj. parameter for distance between 2 forest dislocations for each slip system - atomicVolume, & - t_tr, & !< martensite lamellar thickness [m] for each trans system and instance - p, & !< p-exponent in glide velocity - q, & !< q-exponent in glide velocity - r, & !< r-exponent in twin nucleation rate - s, & !< s-exponent in trans nucleation rate - gamma_char, & !< characteristic shear for twins - B !< drag coefficient - real(pReal), dimension(:,:), allocatable :: & - h_sl_sl, & !< - h_sl_tw, & !< - h_tw_tw, & !< - h_sl_tr, & !< - h_tr_tr !< - integer, dimension(:,:), allocatable :: & - fcc_twinNucleationSlipPair ! ToDo: Better name? Is also use for trans - real(pReal), dimension(:,:), allocatable :: & - n0_sl, & !< slip system normal - forestProjection, & - C66 - real(pReal), dimension(:,:,:), allocatable :: & - P_tr, & - P_sl, & - P_tw, & - C66_tw, & - C66_tr - integer :: & - sum_N_sl, & !< total number of active slip system - sum_N_tw, & !< total number of active twin system - sum_N_tr !< total number of active transformation system - integer, dimension(:), allocatable :: & - N_sl, & !< number of active slip systems for each family - N_tw, & !< number of active twin systems for each family - N_tr !< number of active transformation systems for each family - integer(kind(undefined_ID)), dimension(:), allocatable :: & - outputID !< ID of each post result output - logical :: & - ExtendedDislocations, & !< consider split into partials for climb calculation - fccTwinTransNucleation, & !< twinning and transformation models are for fcc - dipoleFormation !< flag indicating consideration of dipole formation - end type !< container type for internal constitutive parameters - - type :: tDislotwinState - real(pReal), dimension(:,:), pointer :: & - rho_mob, & - rho_dip, & - gamma_sl, & - f_tw, & - f_tr - end type tDislotwinState - - type :: tDislotwinMicrostructure - real(pReal), dimension(:,:), allocatable :: & - Lambda_sl, & !< mean free path between 2 obstacles seen by a moving dislocation - Lambda_tw, & !< mean free path between 2 obstacles seen by a growing twin - Lambda_tr, & !< mean free path between 2 obstacles seen by a growing martensite - tau_pass, & - tau_hat_tw, & - tau_hat_tr, & - V_tw, & !< volume of a new twin - V_tr, & !< volume of a new martensite disc - tau_r_tw, & !< stress to bring partials close together (twin) - tau_r_tr !< stress to bring partials close together (trans) - end type tDislotwinMicrostructure - + + real(pReal), parameter :: & + kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin + + enum, bind(c) + enumerator :: & + undefined_ID, & + rho_mob_ID, & + rho_dip_ID, & + dot_gamma_sl_ID, & + gamma_sl_ID, & + Lambda_sl_ID, & + resolved_stress_slip_ID, & + tau_pass_ID, & + edge_dipole_distance_ID, & + f_tw_ID, & + Lambda_tw_ID, & + resolved_stress_twin_ID, & + tau_hat_tw_ID, & + f_tr_ID + end enum + + type :: tParameters + real(pReal) :: & + mu, & + nu, & + D0, & !< prefactor for self-diffusion coefficient + Qsd, & !< activation energy for dislocation climb + omega, & !< frequency factor for dislocation climb + D, & !< grain size + p_sb, & !< p-exponent in shear band velocity + q_sb, & !< q-exponent in shear band velocity + CEdgeDipMinDistance, & !< + i_tw, & !< + tau_0, & !< strength due to elements in solid solution + L_tw, & !< Length of twin nuclei in Burgers vectors + L_tr, & !< Length of trans nuclei in Burgers vectors + xc_twin, & !< critical distance for formation of twin nucleus + xc_trans, & !< critical distance for formation of trans nucleus + V_cs, & !< cross slip volume + sbResistance, & !< value for shearband resistance (might become an internal state variable at some point) + sbVelocity, & !< value for shearband velocity_0 + sbQedge, & !< activation energy for shear bands + SFE_0K, & !< stacking fault energy at zero K + dSFE_dT, & !< temperature dependance of stacking fault energy + aTol_rho, & !< absolute tolerance for integration of dislocation density + aTol_f_tw, & !< absolute tolerance for integration of twin volume fraction + aTol_f_tr, & !< absolute tolerance for integration of trans volume fraction + gamma_fcc_hex, & !< Free energy difference between austensite and martensite + i_tr, & !< + h !< Stack height of hex nucleus + real(pReal), dimension(:), allocatable :: & + rho_mob_0, & !< initial unipolar dislocation density per slip system + rho_dip_0, & !< initial dipole dislocation density per slip system + b_sl, & !< absolute length of burgers vector [m] for each slip system + b_tw, & !< absolute length of burgers vector [m] for each twin system + b_tr, & !< absolute length of burgers vector [m] for each transformation system + Delta_F,& !< activation energy for glide [J] for each slip system + v0, & !< dislocation velocity prefactor [m/s] for each slip system + dot_N_0_tw, & !< twin nucleation rate [1/m³s] for each twin system + dot_N_0_tr, & !< trans nucleation rate [1/m³s] for each trans system + t_tw, & !< twin thickness [m] for each twin system + CLambdaSlip, & !< Adj. parameter for distance between 2 forest dislocations for each slip system + atomicVolume, & + t_tr, & !< martensite lamellar thickness [m] for each trans system and instance + p, & !< p-exponent in glide velocity + q, & !< q-exponent in glide velocity + r, & !< r-exponent in twin nucleation rate + s, & !< s-exponent in trans nucleation rate + gamma_char, & !< characteristic shear for twins + B !< drag coefficient + real(pReal), dimension(:,:), allocatable :: & + h_sl_sl, & !< + h_sl_tw, & !< + h_tw_tw, & !< + h_sl_tr, & !< + h_tr_tr !< + integer, dimension(:,:), allocatable :: & + fcc_twinNucleationSlipPair ! ToDo: Better name? Is also use for trans + real(pReal), dimension(:,:), allocatable :: & + n0_sl, & !< slip system normal + forestProjection, & + C66 + real(pReal), dimension(:,:,:), allocatable :: & + P_tr, & + P_sl, & + P_tw, & + C66_tw, & + C66_tr + integer :: & + sum_N_sl, & !< total number of active slip system + sum_N_tw, & !< total number of active twin system + sum_N_tr !< total number of active transformation system + integer, dimension(:), allocatable :: & + N_sl, & !< number of active slip systems for each family + N_tw, & !< number of active twin systems for each family + N_tr !< number of active transformation systems for each family + integer(kind(undefined_ID)), dimension(:), allocatable :: & + outputID !< ID of each post result output + logical :: & + ExtendedDislocations, & !< consider split into partials for climb calculation + fccTwinTransNucleation, & !< twinning and transformation models are for fcc + dipoleFormation !< flag indicating consideration of dipole formation + end type !< container type for internal constitutive parameters + + type :: tDislotwinState + real(pReal), dimension(:,:), pointer :: & + rho_mob, & + rho_dip, & + gamma_sl, & + f_tw, & + f_tr + end type tDislotwinState + + type :: tDislotwinMicrostructure + real(pReal), dimension(:,:), allocatable :: & + Lambda_sl, & !< mean free path between 2 obstacles seen by a moving dislocation + Lambda_tw, & !< mean free path between 2 obstacles seen by a growing twin + Lambda_tr, & !< mean free path between 2 obstacles seen by a growing martensite + tau_pass, & + tau_hat_tw, & + tau_hat_tr, & + V_tw, & !< volume of a new twin + V_tr, & !< volume of a new martensite disc + tau_r_tw, & !< stress to bring partials close together (twin) + tau_r_tr !< stress to bring partials close together (trans) + end type tDislotwinMicrostructure + !-------------------------------------------------------------------------------------------------- ! containers for parameters and state - type(tParameters), allocatable, dimension(:) :: param - type(tDislotwinState), allocatable, dimension(:) :: & - dotState, & - state - type(tDislotwinMicrostructure), allocatable, dimension(:) :: dependentState + type(tParameters), allocatable, dimension(:) :: param + type(tDislotwinState), allocatable, dimension(:) :: & + dotState, & + state + type(tDislotwinMicrostructure), allocatable, dimension(:) :: dependentState contains @@ -153,397 +153,397 @@ contains !-------------------------------------------------------------------------------------------------- module subroutine plastic_dislotwin_init - integer :: & - Ninstance, & - p, i, & - NipcMyPhase, outputSize, & - sizeState, sizeDotState, & - startIndex, endIndex - - integer(kind(undefined_ID)) :: & - outputID - - character(len=pStringLen) :: & - extmsg = '' - character(len=pStringLen), dimension(:), allocatable :: & - outputs - - write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_DISLOTWIN_label//' init -+>>>' - - write(6,'(/,a)') ' Ma and Roters, Acta Materialia 52(12):3603–3612, 2004' - write(6,'(a)') ' https://doi.org/10.1016/j.actamat.2004.04.012' - - write(6,'(/,a)') ' Roters et al., Computational Materials Science 39:91–95, 2007' - write(6,'(a)') ' https://doi.org/10.1016/j.commatsci.2006.04.014' - - write(6,'(/,a)') ' Wong et al., Acta Materialia 118:140–151, 2016' - write(6,'(a,/)') ' https://doi.org/10.1016/j.actamat.2016.07.032' - - Ninstance = count(phase_plasticity == PLASTICITY_DISLOTWIN_ID) - - if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & - write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - - allocate(param(Ninstance)) - allocate(state(Ninstance)) - allocate(dotState(Ninstance)) - allocate(dependentState(Ninstance)) - - do p = 1, size(phase_plasticity) - if (phase_plasticity(p) /= PLASTICITY_DISLOTWIN_ID) cycle - associate(prm => param(phase_plasticityInstance(p)), & - dot => dotState(phase_plasticityInstance(p)), & - stt => state(phase_plasticityInstance(p)), & - dst => dependentState(phase_plasticityInstance(p)), & - config => config_phase(p)) - - prm%aTol_rho = config%getFloat('atol_rho', defaultVal=0.0_pReal) - prm%aTol_f_tw = config%getFloat('atol_twinfrac', defaultVal=0.0_pReal) - prm%aTol_f_tr = config%getFloat('atol_transfrac', defaultVal=0.0_pReal) - - ! This data is read in already in lattice - prm%mu = lattice_mu(p) - prm%nu = lattice_nu(p) - prm%C66 = lattice_C66(1:6,1:6,p) - - + integer :: & + Ninstance, & + p, i, & + NipcMyPhase, outputSize, & + sizeState, sizeDotState, & + startIndex, endIndex + + integer(kind(undefined_ID)) :: & + outputID + + character(len=pStringLen) :: & + extmsg = '' + character(len=pStringLen), dimension(:), allocatable :: & + outputs + + write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_DISLOTWIN_label//' init -+>>>' + + write(6,'(/,a)') ' Ma and Roters, Acta Materialia 52(12):3603–3612, 2004' + write(6,'(a)') ' https://doi.org/10.1016/j.actamat.2004.04.012' + + write(6,'(/,a)') ' Roters et al., Computational Materials Science 39:91–95, 2007' + write(6,'(a)') ' https://doi.org/10.1016/j.commatsci.2006.04.014' + + write(6,'(/,a)') ' Wong et al., Acta Materialia 118:140–151, 2016' + write(6,'(a,/)') ' https://doi.org/10.1016/j.actamat.2016.07.032' + + Ninstance = count(phase_plasticity == PLASTICITY_DISLOTWIN_ID) + + if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & + write(6,'(a16,1x,i5,/)') '# instances:',Ninstance + + allocate(param(Ninstance)) + allocate(state(Ninstance)) + allocate(dotState(Ninstance)) + allocate(dependentState(Ninstance)) + + do p = 1, size(phase_plasticity) + if (phase_plasticity(p) /= PLASTICITY_DISLOTWIN_ID) cycle + associate(prm => param(phase_plasticityInstance(p)), & + dot => dotState(phase_plasticityInstance(p)), & + stt => state(phase_plasticityInstance(p)), & + dst => dependentState(phase_plasticityInstance(p)), & + config => config_phase(p)) + + prm%aTol_rho = config%getFloat('atol_rho', defaultVal=0.0_pReal) + prm%aTol_f_tw = config%getFloat('atol_twinfrac', defaultVal=0.0_pReal) + prm%aTol_f_tr = config%getFloat('atol_transfrac', defaultVal=0.0_pReal) + + ! This data is read in already in lattice + prm%mu = lattice_mu(p) + prm%nu = lattice_nu(p) + prm%C66 = lattice_C66(1:6,1:6,p) + + !-------------------------------------------------------------------------------------------------- ! slip related parameters - prm%N_sl = config%getInts('nslip',defaultVal=emptyIntArray) - prm%sum_N_sl = sum(prm%N_sl) - slipActive: if (prm%sum_N_sl > 0) then - prm%P_sl = lattice_SchmidMatrix_slip(prm%N_sl,config%getString('lattice_structure'),& - config%getFloat('c/a',defaultVal=0.0_pReal)) - prm%h_sl_sl = lattice_interaction_SlipBySlip(prm%N_sl, & - config%getFloats('interaction_slipslip'), & - config%getString('lattice_structure')) - prm%forestProjection = lattice_forestProjection_edge(prm%N_sl,config%getString('lattice_structure'),& - config%getFloat('c/a',defaultVal=0.0_pReal)) - prm%forestProjection = transpose(prm%forestProjection) - - prm%n0_sl = lattice_slip_normal(prm%N_sl,config%getString('lattice_structure'),& - config%getFloat('c/a',defaultVal=0.0_pReal)) - prm%fccTwinTransNucleation = merge(.true., .false., lattice_structure(p) == LATTICE_FCC_ID) & - .and. (prm%N_sl(1) == 12) - if(prm%fccTwinTransNucleation) & - prm%fcc_twinNucleationSlipPair = lattice_fcc_twinNucleationSlipPair - - prm%rho_mob_0 = config%getFloats('rhoedge0', requiredSize=size(prm%N_sl)) - prm%rho_dip_0 = config%getFloats('rhoedgedip0',requiredSize=size(prm%N_sl)) - prm%v0 = config%getFloats('v0', requiredSize=size(prm%N_sl)) - prm%b_sl = config%getFloats('slipburgers',requiredSize=size(prm%N_sl)) - prm%Delta_F = config%getFloats('qedge', requiredSize=size(prm%N_sl)) - prm%CLambdaSlip = config%getFloats('clambdaslip',requiredSize=size(prm%N_sl)) - prm%p = config%getFloats('p_slip', requiredSize=size(prm%N_sl)) - prm%q = config%getFloats('q_slip', requiredSize=size(prm%N_sl)) - prm%B = config%getFloats('b', requiredSize=size(prm%N_sl), & - defaultVal=[(0.0_pReal, i=1,size(prm%N_sl))]) - - prm%tau_0 = config%getFloat('solidsolutionstrength') - prm%CEdgeDipMinDistance = config%getFloat('cedgedipmindistance') - prm%D0 = config%getFloat('d0') - prm%Qsd = config%getFloat('qsd') - prm%atomicVolume = config%getFloat('catomicvolume') * prm%b_sl**3.0_pReal - prm%ExtendedDislocations = config%keyExists('/extend_dislocations/') - if (prm%ExtendedDislocations) then - prm%SFE_0K = config%getFloat('sfe_0k') - prm%dSFE_dT = config%getFloat('dsfe_dt') - endif - - ! multiplication factor according to crystal structure (nearest neighbors bcc vs fcc/hex) - !@details: Refer: Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981 - prm%omega = config%getFloat('omega', defaultVal = 1000.0_pReal) & - * merge(12.0_pReal, & - 8.0_pReal, & - lattice_structure(p) == LATTICE_FCC_ID .or. lattice_structure(p) == LATTICE_HEX_ID) - - - ! expand: family => system - prm%rho_mob_0 = math_expand(prm%rho_mob_0, prm%N_sl) - prm%rho_dip_0 = math_expand(prm%rho_dip_0, prm%N_sl) - prm%v0 = math_expand(prm%v0, prm%N_sl) - prm%b_sl = math_expand(prm%b_sl, prm%N_sl) - prm%Delta_F = math_expand(prm%Delta_F, prm%N_sl) - prm%CLambdaSlip = math_expand(prm%CLambdaSlip, prm%N_sl) - prm%p = math_expand(prm%p, prm%N_sl) - prm%q = math_expand(prm%q, prm%N_sl) - prm%B = math_expand(prm%B, prm%N_sl) - prm%atomicVolume = math_expand(prm%atomicVolume,prm%N_sl) - - ! sanity checks - if ( prm%D0 <= 0.0_pReal) extmsg = trim(extmsg)//' D0' - if ( prm%Qsd <= 0.0_pReal) extmsg = trim(extmsg)//' Qsd' - if (any(prm%rho_mob_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_mob_0' - if (any(prm%rho_dip_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_dip_0' - if (any(prm%v0 < 0.0_pReal)) extmsg = trim(extmsg)//' v0' - if (any(prm%b_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' b_sl' - if (any(prm%Delta_F <= 0.0_pReal)) extmsg = trim(extmsg)//' Delta_F' - if (any(prm%CLambdaSlip <= 0.0_pReal)) extmsg = trim(extmsg)//' CLambdaSlip' - if (any(prm%B < 0.0_pReal)) extmsg = trim(extmsg)//' B' - if (any(prm%p<=0.0_pReal .or. prm%p>1.0_pReal)) extmsg = trim(extmsg)//' p' - if (any(prm%q< 1.0_pReal .or. prm%q>2.0_pReal)) extmsg = trim(extmsg)//' q' - - else slipActive - allocate(prm%b_sl(0)) - endif slipActive - + prm%N_sl = config%getInts('nslip',defaultVal=emptyIntArray) + prm%sum_N_sl = sum(prm%N_sl) + slipActive: if (prm%sum_N_sl > 0) then + prm%P_sl = lattice_SchmidMatrix_slip(prm%N_sl,config%getString('lattice_structure'),& + config%getFloat('c/a',defaultVal=0.0_pReal)) + prm%h_sl_sl = lattice_interaction_SlipBySlip(prm%N_sl, & + config%getFloats('interaction_slipslip'), & + config%getString('lattice_structure')) + prm%forestProjection = lattice_forestProjection_edge(prm%N_sl,config%getString('lattice_structure'),& + config%getFloat('c/a',defaultVal=0.0_pReal)) + prm%forestProjection = transpose(prm%forestProjection) + + prm%n0_sl = lattice_slip_normal(prm%N_sl,config%getString('lattice_structure'),& + config%getFloat('c/a',defaultVal=0.0_pReal)) + prm%fccTwinTransNucleation = merge(.true., .false., lattice_structure(p) == LATTICE_FCC_ID) & + .and. (prm%N_sl(1) == 12) + if(prm%fccTwinTransNucleation) & + prm%fcc_twinNucleationSlipPair = lattice_fcc_twinNucleationSlipPair + + prm%rho_mob_0 = config%getFloats('rhoedge0', requiredSize=size(prm%N_sl)) + prm%rho_dip_0 = config%getFloats('rhoedgedip0',requiredSize=size(prm%N_sl)) + prm%v0 = config%getFloats('v0', requiredSize=size(prm%N_sl)) + prm%b_sl = config%getFloats('slipburgers',requiredSize=size(prm%N_sl)) + prm%Delta_F = config%getFloats('qedge', requiredSize=size(prm%N_sl)) + prm%CLambdaSlip = config%getFloats('clambdaslip',requiredSize=size(prm%N_sl)) + prm%p = config%getFloats('p_slip', requiredSize=size(prm%N_sl)) + prm%q = config%getFloats('q_slip', requiredSize=size(prm%N_sl)) + prm%B = config%getFloats('b', requiredSize=size(prm%N_sl), & + defaultVal=[(0.0_pReal, i=1,size(prm%N_sl))]) + + prm%tau_0 = config%getFloat('solidsolutionstrength') + prm%CEdgeDipMinDistance = config%getFloat('cedgedipmindistance') + prm%D0 = config%getFloat('d0') + prm%Qsd = config%getFloat('qsd') + prm%atomicVolume = config%getFloat('catomicvolume') * prm%b_sl**3.0_pReal + prm%ExtendedDislocations = config%keyExists('/extend_dislocations/') + if (prm%ExtendedDislocations) then + prm%SFE_0K = config%getFloat('sfe_0k') + prm%dSFE_dT = config%getFloat('dsfe_dt') + endif + + ! multiplication factor according to crystal structure (nearest neighbors bcc vs fcc/hex) + !@details: Refer: Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981 + prm%omega = config%getFloat('omega', defaultVal = 1000.0_pReal) & + * merge(12.0_pReal, & + 8.0_pReal, & + lattice_structure(p) == LATTICE_FCC_ID .or. lattice_structure(p) == LATTICE_HEX_ID) + + + ! expand: family => system + prm%rho_mob_0 = math_expand(prm%rho_mob_0, prm%N_sl) + prm%rho_dip_0 = math_expand(prm%rho_dip_0, prm%N_sl) + prm%v0 = math_expand(prm%v0, prm%N_sl) + prm%b_sl = math_expand(prm%b_sl, prm%N_sl) + prm%Delta_F = math_expand(prm%Delta_F, prm%N_sl) + prm%CLambdaSlip = math_expand(prm%CLambdaSlip, prm%N_sl) + prm%p = math_expand(prm%p, prm%N_sl) + prm%q = math_expand(prm%q, prm%N_sl) + prm%B = math_expand(prm%B, prm%N_sl) + prm%atomicVolume = math_expand(prm%atomicVolume,prm%N_sl) + + ! sanity checks + if ( prm%D0 <= 0.0_pReal) extmsg = trim(extmsg)//' D0' + if ( prm%Qsd <= 0.0_pReal) extmsg = trim(extmsg)//' Qsd' + if (any(prm%rho_mob_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_mob_0' + if (any(prm%rho_dip_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_dip_0' + if (any(prm%v0 < 0.0_pReal)) extmsg = trim(extmsg)//' v0' + if (any(prm%b_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' b_sl' + if (any(prm%Delta_F <= 0.0_pReal)) extmsg = trim(extmsg)//' Delta_F' + if (any(prm%CLambdaSlip <= 0.0_pReal)) extmsg = trim(extmsg)//' CLambdaSlip' + if (any(prm%B < 0.0_pReal)) extmsg = trim(extmsg)//' B' + if (any(prm%p<=0.0_pReal .or. prm%p>1.0_pReal)) extmsg = trim(extmsg)//' p' + if (any(prm%q< 1.0_pReal .or. prm%q>2.0_pReal)) extmsg = trim(extmsg)//' q' + + else slipActive + allocate(prm%b_sl(0)) + endif slipActive + !-------------------------------------------------------------------------------------------------- ! twin related parameters - prm%N_tw = config%getInts('ntwin', defaultVal=emptyIntArray) - prm%sum_N_tw = sum(prm%N_tw) - if (prm%sum_N_tw > 0) then - prm%P_tw = lattice_SchmidMatrix_twin(prm%N_tw,config%getString('lattice_structure'),& - config%getFloat('c/a',defaultVal=0.0_pReal)) - prm%h_tw_tw = lattice_interaction_TwinByTwin(prm%N_tw,& - config%getFloats('interaction_twintwin'), & - config%getString('lattice_structure')) - - prm%b_tw = config%getFloats('twinburgers', requiredSize=size(prm%N_tw)) - prm%t_tw = config%getFloats('twinsize', requiredSize=size(prm%N_tw)) - prm%r = config%getFloats('r_twin', requiredSize=size(prm%N_tw)) - - prm%xc_twin = config%getFloat('xc_twin') - prm%L_tw = config%getFloat('l0_twin') - prm%i_tw = config%getFloat('cmfptwin') - - prm%gamma_char= lattice_characteristicShear_Twin(prm%N_tw,config%getString('lattice_structure'),& - config%getFloat('c/a',defaultVal=0.0_pReal)) - - prm%C66_tw = lattice_C66_twin(prm%N_tw,prm%C66,config%getString('lattice_structure'),& - config%getFloat('c/a',defaultVal=0.0_pReal)) - - if (.not. prm%fccTwinTransNucleation) then - prm%dot_N_0_tw = config%getFloats('ndot0_twin') - prm%dot_N_0_tw = math_expand(prm%dot_N_0_tw,prm%N_tw) - endif - - ! expand: family => system - prm%b_tw = math_expand(prm%b_tw,prm%N_tw) - prm%t_tw = math_expand(prm%t_tw,prm%N_tw) - prm%r = math_expand(prm%r,prm%N_tw) - - else - allocate(prm%gamma_char(0)) - allocate(prm%t_tw (0)) - allocate(prm%b_tw (0)) - allocate(prm%r (0)) - allocate(prm%h_tw_tw (0,0)) - endif - + prm%N_tw = config%getInts('ntwin', defaultVal=emptyIntArray) + prm%sum_N_tw = sum(prm%N_tw) + if (prm%sum_N_tw > 0) then + prm%P_tw = lattice_SchmidMatrix_twin(prm%N_tw,config%getString('lattice_structure'),& + config%getFloat('c/a',defaultVal=0.0_pReal)) + prm%h_tw_tw = lattice_interaction_TwinByTwin(prm%N_tw,& + config%getFloats('interaction_twintwin'), & + config%getString('lattice_structure')) + + prm%b_tw = config%getFloats('twinburgers', requiredSize=size(prm%N_tw)) + prm%t_tw = config%getFloats('twinsize', requiredSize=size(prm%N_tw)) + prm%r = config%getFloats('r_twin', requiredSize=size(prm%N_tw)) + + prm%xc_twin = config%getFloat('xc_twin') + prm%L_tw = config%getFloat('l0_twin') + prm%i_tw = config%getFloat('cmfptwin') + + prm%gamma_char= lattice_characteristicShear_Twin(prm%N_tw,config%getString('lattice_structure'),& + config%getFloat('c/a',defaultVal=0.0_pReal)) + + prm%C66_tw = lattice_C66_twin(prm%N_tw,prm%C66,config%getString('lattice_structure'),& + config%getFloat('c/a',defaultVal=0.0_pReal)) + + if (.not. prm%fccTwinTransNucleation) then + prm%dot_N_0_tw = config%getFloats('ndot0_twin') + prm%dot_N_0_tw = math_expand(prm%dot_N_0_tw,prm%N_tw) + endif + + ! expand: family => system + prm%b_tw = math_expand(prm%b_tw,prm%N_tw) + prm%t_tw = math_expand(prm%t_tw,prm%N_tw) + prm%r = math_expand(prm%r,prm%N_tw) + + else + allocate(prm%gamma_char(0)) + allocate(prm%t_tw (0)) + allocate(prm%b_tw (0)) + allocate(prm%r (0)) + allocate(prm%h_tw_tw (0,0)) + endif + !-------------------------------------------------------------------------------------------------- ! transformation related parameters - prm%N_tr = config%getInts('ntrans', defaultVal=emptyIntArray) - prm%sum_N_tr = sum(prm%N_tr) - if (prm%sum_N_tr > 0) then - prm%b_tr = config%getFloats('transburgers') - prm%b_tr = math_expand(prm%b_tr,prm%N_tr) - - prm%h = config%getFloat('transstackheight', defaultVal=0.0_pReal) ! ToDo: How to handle that??? - prm%i_tr = config%getFloat('cmfptrans', defaultVal=0.0_pReal) ! ToDo: How to handle that??? - prm%gamma_fcc_hex = config%getFloat('deltag') - prm%xc_trans = config%getFloat('xc_trans', defaultVal=0.0_pReal) ! ToDo: How to handle that??? - prm%L_tr = config%getFloat('l0_trans') - - prm%h_tr_tr = lattice_interaction_TransByTrans(prm%N_tr,& - config%getFloats('interaction_transtrans'), & - config%getString('lattice_structure')) - - prm%C66_tr = lattice_C66_trans(prm%N_tr,prm%C66, & - config%getString('trans_lattice_structure'), & - 0.0_pReal, & - config%getFloat('a_bcc', defaultVal=0.0_pReal), & - config%getFloat('a_fcc', defaultVal=0.0_pReal)) - - prm%P_tr = lattice_SchmidMatrix_trans(prm%N_tr, & - config%getString('trans_lattice_structure'), & - 0.0_pReal, & - config%getFloat('a_bcc', defaultVal=0.0_pReal), & - config%getFloat('a_fcc', defaultVal=0.0_pReal)) - - if (lattice_structure(p) /= LATTICE_fcc_ID) then - prm%dot_N_0_tr = config%getFloats('ndot0_trans') - prm%dot_N_0_tr = math_expand(prm%dot_N_0_tr,prm%N_tr) - endif - prm%t_tr = config%getFloats('lamellarsize') - prm%t_tr = math_expand(prm%t_tr,prm%N_tr) - prm%s = config%getFloats('s_trans',defaultVal=[0.0_pReal]) - prm%s = math_expand(prm%s,prm%N_tr) - else - allocate(prm%t_tr (0)) - allocate(prm%b_tr (0)) - allocate(prm%s (0)) - allocate(prm%h_tr_tr(0,0)) - endif + prm%N_tr = config%getInts('ntrans', defaultVal=emptyIntArray) + prm%sum_N_tr = sum(prm%N_tr) + if (prm%sum_N_tr > 0) then + prm%b_tr = config%getFloats('transburgers') + prm%b_tr = math_expand(prm%b_tr,prm%N_tr) + + prm%h = config%getFloat('transstackheight', defaultVal=0.0_pReal) ! ToDo: How to handle that??? + prm%i_tr = config%getFloat('cmfptrans', defaultVal=0.0_pReal) ! ToDo: How to handle that??? + prm%gamma_fcc_hex = config%getFloat('deltag') + prm%xc_trans = config%getFloat('xc_trans', defaultVal=0.0_pReal) ! ToDo: How to handle that??? + prm%L_tr = config%getFloat('l0_trans') + + prm%h_tr_tr = lattice_interaction_TransByTrans(prm%N_tr,& + config%getFloats('interaction_transtrans'), & + config%getString('lattice_structure')) + + prm%C66_tr = lattice_C66_trans(prm%N_tr,prm%C66, & + config%getString('trans_lattice_structure'), & + 0.0_pReal, & + config%getFloat('a_bcc', defaultVal=0.0_pReal), & + config%getFloat('a_fcc', defaultVal=0.0_pReal)) + + prm%P_tr = lattice_SchmidMatrix_trans(prm%N_tr, & + config%getString('trans_lattice_structure'), & + 0.0_pReal, & + config%getFloat('a_bcc', defaultVal=0.0_pReal), & + config%getFloat('a_fcc', defaultVal=0.0_pReal)) + + if (lattice_structure(p) /= LATTICE_fcc_ID) then + prm%dot_N_0_tr = config%getFloats('ndot0_trans') + prm%dot_N_0_tr = math_expand(prm%dot_N_0_tr,prm%N_tr) + endif + prm%t_tr = config%getFloats('lamellarsize') + prm%t_tr = math_expand(prm%t_tr,prm%N_tr) + prm%s = config%getFloats('s_trans',defaultVal=[0.0_pReal]) + prm%s = math_expand(prm%s,prm%N_tr) + else + allocate(prm%t_tr (0)) + allocate(prm%b_tr (0)) + allocate(prm%s (0)) + allocate(prm%h_tr_tr(0,0)) + endif + + if (sum(prm%N_tw) > 0 .or. prm%sum_N_tr > 0) then + prm%SFE_0K = config%getFloat('sfe_0k') + prm%dSFE_dT = config%getFloat('dsfe_dt') + prm%V_cs = config%getFloat('vcrossslip') + endif + + if (prm%sum_N_sl > 0 .and. prm%sum_N_tw > 0) then + prm%h_sl_tw = lattice_interaction_SlipByTwin(prm%N_sl,prm%N_tw,& + config%getFloats('interaction_sliptwin'), & + config%getString('lattice_structure')) + if (prm%fccTwinTransNucleation .and. prm%sum_N_tw > 12) write(6,*) 'mist' ! ToDo: implement better test. The model will fail also if N_tw is [6,6] + endif + + if (prm%sum_N_sl > 0 .and. prm%sum_N_tr > 0) then + prm%h_sl_tr = lattice_interaction_SlipByTrans(prm%N_sl,prm%N_tr,& + config%getFloats('interaction_sliptrans'), & + config%getString('lattice_structure')) + if (prm%fccTwinTransNucleation .and. prm%sum_N_tr > 12) write(6,*) 'mist' ! ToDo: implement better test. The model will fail also if N_tr is [6,6] + endif - if (sum(prm%N_tw) > 0 .or. prm%sum_N_tr > 0) then - prm%SFE_0K = config%getFloat('sfe_0k') - prm%dSFE_dT = config%getFloat('dsfe_dt') - prm%V_cs = config%getFloat('vcrossslip') - endif - - if (prm%sum_N_sl > 0 .and. prm%sum_N_tw > 0) then - prm%h_sl_tw = lattice_interaction_SlipByTwin(prm%N_sl,prm%N_tw,& - config%getFloats('interaction_sliptwin'), & - config%getString('lattice_structure')) - if (prm%fccTwinTransNucleation .and. prm%sum_N_tw > 12) write(6,*) 'mist' ! ToDo: implement better test. The model will fail also if N_tw is [6,6] - endif - - if (prm%sum_N_sl > 0 .and. prm%sum_N_tr > 0) then - prm%h_sl_tr = lattice_interaction_SlipByTrans(prm%N_sl,prm%N_tr,& - config%getFloats('interaction_sliptrans'), & - config%getString('lattice_structure')) - if (prm%fccTwinTransNucleation .and. prm%sum_N_tr > 12) write(6,*) 'mist' ! ToDo: implement better test. The model will fail also if N_tr is [6,6] - endif - !-------------------------------------------------------------------------------------------------- ! shearband related parameters - prm%sbVelocity = config%getFloat('shearbandvelocity',defaultVal=0.0_pReal) - if (prm%sbVelocity > 0.0_pReal) then - prm%sbResistance = config%getFloat('shearbandresistance') - prm%sbQedge = config%getFloat('qedgepersbsystem') - prm%p_sb = config%getFloat('p_shearband') - prm%q_sb = config%getFloat('q_shearband') - - ! sanity checks - if (prm%sbResistance < 0.0_pReal) extmsg = trim(extmsg)//' shearbandresistance' - if (prm%sbQedge < 0.0_pReal) extmsg = trim(extmsg)//' qedgepersbsystem' - if (prm%p_sb <= 0.0_pReal) extmsg = trim(extmsg)//' p_shearband' - if (prm%q_sb <= 0.0_pReal) extmsg = trim(extmsg)//' q_shearband' - endif - - - - prm%D = config%getFloat('grainsize') - - if (config%keyExists('dipoleformationfactor')) call IO_error(1,ext_msg='use /nodipoleformation/') - prm%dipoleformation = .not. config%keyExists('/nodipoleformation/') - - - !if (Ndot0PerTwinFamily(f,p) < 0.0_pReal) & - ! call IO_error(211,el=p,ext_msg='dot_N_0_tw ('//PLASTICITY_DISLOTWIN_label//')') - - if (any(prm%atomicVolume <= 0.0_pReal)) & - call IO_error(211,el=p,ext_msg='cAtomicVolume ('//PLASTICITY_DISLOTWIN_label//')') - if (prm%sum_N_tw > 0) then - if (prm%aTol_rho <= 0.0_pReal) & - call IO_error(211,el=p,ext_msg='aTol_rho ('//PLASTICITY_DISLOTWIN_label//')') - if (prm%aTol_f_tw <= 0.0_pReal) & - call IO_error(211,el=p,ext_msg='aTol_f_tw ('//PLASTICITY_DISLOTWIN_label//')') - endif - if (prm%sum_N_tr > 0) then - if (prm%aTol_f_tr <= 0.0_pReal) & - call IO_error(211,el=p,ext_msg='aTol_f_tr ('//PLASTICITY_DISLOTWIN_label//')') - endif + prm%sbVelocity = config%getFloat('shearbandvelocity',defaultVal=0.0_pReal) + if (prm%sbVelocity > 0.0_pReal) then + prm%sbResistance = config%getFloat('shearbandresistance') + prm%sbQedge = config%getFloat('qedgepersbsystem') + prm%p_sb = config%getFloat('p_shearband') + prm%q_sb = config%getFloat('q_shearband') + + ! sanity checks + if (prm%sbResistance < 0.0_pReal) extmsg = trim(extmsg)//' shearbandresistance' + if (prm%sbQedge < 0.0_pReal) extmsg = trim(extmsg)//' qedgepersbsystem' + if (prm%p_sb <= 0.0_pReal) extmsg = trim(extmsg)//' p_shearband' + if (prm%q_sb <= 0.0_pReal) extmsg = trim(extmsg)//' q_shearband' + endif - outputs = config%getStrings('(output)', defaultVal=emptyStringArray) - allocate(prm%outputID(0)) - do i= 1, size(outputs) - outputID = undefined_ID - select case(outputs(i)) - case ('rho_mob') - outputID = merge(rho_mob_ID,undefined_ID,prm%sum_N_sl > 0) - outputSize = prm%sum_N_sl - case ('rho_dip') - outputID = merge(rho_dip_ID,undefined_ID,prm%sum_N_sl > 0) - outputSize = prm%sum_N_sl - case ('gamma_sl') - outputID = merge(gamma_sl_ID,undefined_ID,prm%sum_N_sl > 0) - outputSize = prm%sum_N_sl - case ('lambda_sl') - outputID = merge(Lambda_sl_ID,undefined_ID,prm%sum_N_sl > 0) - outputSize = prm%sum_N_sl - case ('tau_pass') - outputID= merge(tau_pass_ID,undefined_ID,prm%sum_N_sl > 0) - outputSize = prm%sum_N_sl - - case ('f_tw') - outputID = merge(f_tw_ID,undefined_ID,prm%sum_N_tw >0) - outputSize = prm%sum_N_tw - case ('lambda_tw') - outputID = merge(Lambda_tw_ID,undefined_ID,prm%sum_N_tw >0) - outputSize = prm%sum_N_tw - case ('tau_hat_tw') - outputID = merge(tau_hat_tw_ID,undefined_ID,prm%sum_N_tw >0) - outputSize = prm%sum_N_tw + + + prm%D = config%getFloat('grainsize') + + if (config%keyExists('dipoleformationfactor')) call IO_error(1,ext_msg='use /nodipoleformation/') + prm%dipoleformation = .not. config%keyExists('/nodipoleformation/') + + + !if (Ndot0PerTwinFamily(f,p) < 0.0_pReal) & + ! call IO_error(211,el=p,ext_msg='dot_N_0_tw ('//PLASTICITY_DISLOTWIN_label//')') + + if (any(prm%atomicVolume <= 0.0_pReal)) & + call IO_error(211,el=p,ext_msg='cAtomicVolume ('//PLASTICITY_DISLOTWIN_label//')') + if (prm%sum_N_tw > 0) then + if (prm%aTol_rho <= 0.0_pReal) & + call IO_error(211,el=p,ext_msg='aTol_rho ('//PLASTICITY_DISLOTWIN_label//')') + if (prm%aTol_f_tw <= 0.0_pReal) & + call IO_error(211,el=p,ext_msg='aTol_f_tw ('//PLASTICITY_DISLOTWIN_label//')') + endif + if (prm%sum_N_tr > 0) then + if (prm%aTol_f_tr <= 0.0_pReal) & + call IO_error(211,el=p,ext_msg='aTol_f_tr ('//PLASTICITY_DISLOTWIN_label//')') + endif + + outputs = config%getStrings('(output)', defaultVal=emptyStringArray) + allocate(prm%outputID(0)) + do i= 1, size(outputs) + outputID = undefined_ID + select case(outputs(i)) + case ('rho_mob') + outputID = merge(rho_mob_ID,undefined_ID,prm%sum_N_sl > 0) + outputSize = prm%sum_N_sl + case ('rho_dip') + outputID = merge(rho_dip_ID,undefined_ID,prm%sum_N_sl > 0) + outputSize = prm%sum_N_sl + case ('gamma_sl') + outputID = merge(gamma_sl_ID,undefined_ID,prm%sum_N_sl > 0) + outputSize = prm%sum_N_sl + case ('lambda_sl') + outputID = merge(Lambda_sl_ID,undefined_ID,prm%sum_N_sl > 0) + outputSize = prm%sum_N_sl + case ('tau_pass') + outputID= merge(tau_pass_ID,undefined_ID,prm%sum_N_sl > 0) + outputSize = prm%sum_N_sl + + case ('f_tw') + outputID = merge(f_tw_ID,undefined_ID,prm%sum_N_tw >0) + outputSize = prm%sum_N_tw + case ('lambda_tw') + outputID = merge(Lambda_tw_ID,undefined_ID,prm%sum_N_tw >0) + outputSize = prm%sum_N_tw + case ('tau_hat_tw') + outputID = merge(tau_hat_tw_ID,undefined_ID,prm%sum_N_tw >0) + outputSize = prm%sum_N_tw + + case ('f_tr') + outputID = f_tr_ID + outputSize = prm%sum_N_tr - case ('f_tr') - outputID = f_tr_ID - outputSize = prm%sum_N_tr - - end select - - if (outputID /= undefined_ID) then - prm%outputID = [prm%outputID, outputID] - endif - - enddo - + end select + + if (outputID /= undefined_ID) then + prm%outputID = [prm%outputID, outputID] + endif + + enddo + !-------------------------------------------------------------------------------------------------- ! allocate state arrays - NipcMyPhase = count(material_phaseAt == p) * discretization_nIP - 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 - sizeState = sizeDotState - - call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0) - + NipcMyPhase = count(material_phaseAt == p) * discretization_nIP + 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 + sizeState = sizeDotState + + call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0) + !-------------------------------------------------------------------------------------------------- ! locally defined state aliases and initialization of state0 and aTolState - startIndex = 1 - endIndex = prm%sum_N_sl - stt%rho_mob=>plasticState(p)%state(startIndex:endIndex,:) - stt%rho_mob= spread(prm%rho_mob_0,2,NipcMyPhase) - dot%rho_mob=>plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%aTolState(startIndex:endIndex) = prm%aTol_rho - - startIndex = endIndex + 1 - endIndex = endIndex + prm%sum_N_sl - stt%rho_dip=>plasticState(p)%state(startIndex:endIndex,:) - stt%rho_dip= spread(prm%rho_dip_0,2,NipcMyPhase) - dot%rho_dip=>plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%aTolState(startIndex:endIndex) = prm%aTol_rho - - startIndex = endIndex + 1 - endIndex = endIndex + prm%sum_N_sl - stt%gamma_sl=>plasticState(p)%state(startIndex:endIndex,:) - dot%gamma_sl=>plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%aTolState(startIndex:endIndex) = 1.0e6_pReal !ToDo: better make optional parameter - ! global alias - plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%accumulatedSlip => plasticState(p)%state(startIndex:endIndex,:) - - startIndex = endIndex + 1 - endIndex = endIndex + prm%sum_N_tw - stt%f_tw=>plasticState(p)%state(startIndex:endIndex,:) - dot%f_tw=>plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%aTolState(startIndex:endIndex) = prm%aTol_f_tw - - startIndex = endIndex + 1 - endIndex = endIndex + prm%sum_N_tr - stt%f_tr=>plasticState(p)%state(startIndex:endIndex,:) - dot%f_tr=>plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%aTolState(startIndex:endIndex) = prm%aTol_f_tr - - allocate(dst%Lambda_sl (prm%sum_N_sl,NipcMyPhase),source=0.0_pReal) - allocate(dst%tau_pass (prm%sum_N_sl,NipcMyPhase),source=0.0_pReal) - - allocate(dst%Lambda_tw (prm%sum_N_tw,NipcMyPhase),source=0.0_pReal) - allocate(dst%tau_hat_tw (prm%sum_N_tw,NipcMyPhase),source=0.0_pReal) - allocate(dst%tau_r_tw (prm%sum_N_tw,NipcMyPhase),source=0.0_pReal) - allocate(dst%V_tw (prm%sum_N_tw,NipcMyPhase),source=0.0_pReal) - - allocate(dst%Lambda_tr (prm%sum_N_tr,NipcMyPhase),source=0.0_pReal) - allocate(dst%tau_hat_tr (prm%sum_N_tr,NipcMyPhase),source=0.0_pReal) - allocate(dst%tau_r_tr (prm%sum_N_tr,NipcMyPhase),source=0.0_pReal) - allocate(dst%V_tr (prm%sum_N_tr,NipcMyPhase),source=0.0_pReal) - - - plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally - - end associate - - enddo + startIndex = 1 + endIndex = prm%sum_N_sl + stt%rho_mob=>plasticState(p)%state(startIndex:endIndex,:) + stt%rho_mob= spread(prm%rho_mob_0,2,NipcMyPhase) + dot%rho_mob=>plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%aTolState(startIndex:endIndex) = prm%aTol_rho + + startIndex = endIndex + 1 + endIndex = endIndex + prm%sum_N_sl + stt%rho_dip=>plasticState(p)%state(startIndex:endIndex,:) + stt%rho_dip= spread(prm%rho_dip_0,2,NipcMyPhase) + dot%rho_dip=>plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%aTolState(startIndex:endIndex) = prm%aTol_rho + + startIndex = endIndex + 1 + endIndex = endIndex + prm%sum_N_sl + stt%gamma_sl=>plasticState(p)%state(startIndex:endIndex,:) + dot%gamma_sl=>plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%aTolState(startIndex:endIndex) = 1.0e6_pReal !ToDo: better make optional parameter + ! global alias + plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%accumulatedSlip => plasticState(p)%state(startIndex:endIndex,:) + + startIndex = endIndex + 1 + endIndex = endIndex + prm%sum_N_tw + stt%f_tw=>plasticState(p)%state(startIndex:endIndex,:) + dot%f_tw=>plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%aTolState(startIndex:endIndex) = prm%aTol_f_tw + + startIndex = endIndex + 1 + endIndex = endIndex + prm%sum_N_tr + stt%f_tr=>plasticState(p)%state(startIndex:endIndex,:) + dot%f_tr=>plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%aTolState(startIndex:endIndex) = prm%aTol_f_tr + + allocate(dst%Lambda_sl (prm%sum_N_sl,NipcMyPhase),source=0.0_pReal) + allocate(dst%tau_pass (prm%sum_N_sl,NipcMyPhase),source=0.0_pReal) + + allocate(dst%Lambda_tw (prm%sum_N_tw,NipcMyPhase),source=0.0_pReal) + allocate(dst%tau_hat_tw (prm%sum_N_tw,NipcMyPhase),source=0.0_pReal) + allocate(dst%tau_r_tw (prm%sum_N_tw,NipcMyPhase),source=0.0_pReal) + allocate(dst%V_tw (prm%sum_N_tw,NipcMyPhase),source=0.0_pReal) + + allocate(dst%Lambda_tr (prm%sum_N_tr,NipcMyPhase),source=0.0_pReal) + allocate(dst%tau_hat_tr (prm%sum_N_tr,NipcMyPhase),source=0.0_pReal) + allocate(dst%tau_r_tr (prm%sum_N_tr,NipcMyPhase),source=0.0_pReal) + allocate(dst%V_tr (prm%sum_N_tr,NipcMyPhase),source=0.0_pReal) + + + plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally + + end associate + + enddo end subroutine plastic_dislotwin_init @@ -553,36 +553,36 @@ end subroutine plastic_dislotwin_init !-------------------------------------------------------------------------------------------------- module function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC) - real(pReal), dimension(6,6) :: & - homogenizedC - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - - integer :: i, & - of - real(pReal) :: f_unrotated - - of = material_phasememberAt(ipc,ip,el) - associate(prm => param(phase_plasticityInstance(material_phaseAt(ipc,el))),& - stt => state(phase_plasticityInstance(material_phaseAT(ipc,el)))) - - 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)) - - homogenizedC = f_unrotated * prm%C66 - do i=1,prm%sum_N_tw - homogenizedC = homogenizedC & - + stt%f_tw(i,of)*prm%C66_tw(1:6,1:6,i) - enddo - do i=1,prm%sum_N_tr - homogenizedC = homogenizedC & - + stt%f_tr(i,of)*prm%C66_tr(1:6,1:6,i) - enddo - - end associate + real(pReal), dimension(6,6) :: & + homogenizedC + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + + integer :: i, & + of + real(pReal) :: f_unrotated + + of = material_phasememberAt(ipc,ip,el) + associate(prm => param(phase_plasticityInstance(material_phaseAt(ipc,el))),& + stt => state(phase_plasticityInstance(material_phaseAT(ipc,el)))) + + 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)) + + homogenizedC = f_unrotated * prm%C66 + do i=1,prm%sum_N_tw + homogenizedC = homogenizedC & + + stt%f_tw(i,of)*prm%C66_tw(1:6,1:6,i) + enddo + do i=1,prm%sum_N_tr + homogenizedC = homogenizedC & + + stt%f_tr(i,of)*prm%C66_tr(1:6,1:6,i) + enddo + + end associate end function plastic_dislotwin_homogenizedC @@ -592,112 +592,113 @@ end function plastic_dislotwin_homogenizedC !-------------------------------------------------------------------------------------------------- module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) - 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 - real(pReal), intent(in) :: T - - integer :: i,k,l,m,n - real(pReal) :: f_unrotated,StressRatio_p,& - BoltzmannRatio, & - ddot_gamma_dtau, & - tau - real(pReal), dimension(param(instance)%sum_N_sl) :: & + 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 + real(pReal), intent(in) :: T + + integer :: i,k,l,m,n + real(pReal) :: & + f_unrotated,StressRatio_p,& + BoltzmannRatio, & + ddot_gamma_dtau, & + tau + real(pReal), dimension(param(instance)%sum_N_sl) :: & dot_gamma_sl,ddot_gamma_dtau_slip - real(pReal), dimension(param(instance)%sum_N_tw) :: & + real(pReal), dimension(param(instance)%sum_N_tw) :: & dot_gamma_twin,ddot_gamma_dtau_twin - real(pReal), dimension(param(instance)%sum_N_tr) :: & + real(pReal), dimension(param(instance)%sum_N_tr) :: & dot_gamma_tr,ddot_gamma_dtau_trans - real(pReal):: dot_gamma_sb - real(pReal), dimension(3,3) :: eigVectors, P_sb - real(pReal), dimension(3) :: eigValues - logical :: error - real(pReal), dimension(3,6), parameter :: & - sb_sComposition = & - reshape(real([& - 1, 0, 1, & - 1, 0,-1, & - 1, 1, 0, & - 1,-1, 0, & - 0, 1, 1, & - 0, 1,-1 & - ],pReal),[ 3,6]), & - sb_mComposition = & - reshape(real([& - 1, 0,-1, & - 1, 0,+1, & - 1,-1, 0, & - 1, 1, 0, & - 0, 1,-1, & - 0, 1, 1 & - ],pReal),[ 3,6]) - - associate(prm => param(instance), stt => state(instance)) - - 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)) - - Lp = 0.0_pReal - dLp_dMp = 0.0_pReal - - call kinetics_slip(Mp,T,instance,of,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) & - dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & - + ddot_gamma_dtau_slip(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i) - enddo slipContribution - - !ToDo: Why do this before shear banding? - Lp = Lp * f_unrotated - dLp_dMp = dLp_dMp * f_unrotated + real(pReal):: dot_gamma_sb + real(pReal), dimension(3,3) :: eigVectors, P_sb + real(pReal), dimension(3) :: eigValues + logical :: error + real(pReal), dimension(3,6), parameter :: & + sb_sComposition = & + reshape(real([& + 1, 0, 1, & + 1, 0,-1, & + 1, 1, 0, & + 1,-1, 0, & + 0, 1, 1, & + 0, 1,-1 & + ],pReal),[ 3,6]), & + sb_mComposition = & + reshape(real([& + 1, 0,-1, & + 1, 0,+1, & + 1,-1, 0, & + 1, 1, 0, & + 0, 1,-1, & + 0, 1, 1 & + ],pReal),[ 3,6]) - shearBandingContribution: if(dNeq0(prm%sbVelocity)) then - - BoltzmannRatio = prm%sbQedge/(kB*T) - call math_eigenValuesVectorsSym(Mp,eigValues,eigVectors,error) - - do i = 1,6 - P_sb = 0.5_pReal * math_outer(matmul(eigVectors,sb_sComposition(1:3,i)),& - matmul(eigVectors,sb_mComposition(1:3,i))) - tau = math_mul33xx33(Mp,P_sb) - - significantShearBandStress: if (abs(tau) > tol_math_check) then - StressRatio_p = (abs(tau)/prm%sbResistance)**prm%p_sb - dot_gamma_sb = sign(prm%sbVelocity*exp(-BoltzmannRatio*(1-StressRatio_p)**prm%q_sb), tau) - ddot_gamma_dtau = abs(dot_gamma_sb)*BoltzmannRatio* prm%p_sb*prm%q_sb/ prm%sbResistance & - * (abs(tau)/prm%sbResistance)**(prm%p_sb-1.0_pReal) & - * (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal) + associate(prm => param(instance), stt => state(instance)) - Lp = Lp + dot_gamma_sb * P_sb - forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & - + ddot_gamma_dtau * P_sb(k,l) * P_sb(m,n) - endif significantShearBandStress - enddo - - endif shearBandingContribution + 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)) - call kinetics_twin(Mp,T,dot_gamma_sl,instance,of,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) * f_unrotated - forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & - + ddot_gamma_dtau_twin(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i) * f_unrotated - enddo twinContibution + Lp = 0.0_pReal + dLp_dMp = 0.0_pReal - call kinetics_trans(Mp,T,dot_gamma_sl,instance,of,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) * f_unrotated - forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & - + ddot_gamma_dtau_trans(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i) * f_unrotated - enddo transContibution - - - end associate + call kinetics_slip(Mp,T,instance,of,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) & + dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + + ddot_gamma_dtau_slip(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i) + enddo slipContribution + + !ToDo: Why do this before shear banding? + Lp = Lp * f_unrotated + dLp_dMp = dLp_dMp * f_unrotated + + shearBandingContribution: if(dNeq0(prm%sbVelocity)) then + + BoltzmannRatio = prm%sbQedge/(kB*T) + call math_eigenValuesVectorsSym(Mp,eigValues,eigVectors,error) + + do i = 1,6 + P_sb = 0.5_pReal * math_outer(matmul(eigVectors,sb_sComposition(1:3,i)),& + matmul(eigVectors,sb_mComposition(1:3,i))) + tau = math_mul33xx33(Mp,P_sb) + + significantShearBandStress: if (abs(tau) > tol_math_check) then + StressRatio_p = (abs(tau)/prm%sbResistance)**prm%p_sb + dot_gamma_sb = sign(prm%sbVelocity*exp(-BoltzmannRatio*(1-StressRatio_p)**prm%q_sb), tau) + ddot_gamma_dtau = abs(dot_gamma_sb)*BoltzmannRatio* prm%p_sb*prm%q_sb/ prm%sbResistance & + * (abs(tau)/prm%sbResistance)**(prm%p_sb-1.0_pReal) & + * (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal) + + Lp = Lp + dot_gamma_sb * P_sb + forall (k=1:3,l=1:3,m=1:3,n=1:3) & + dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + + ddot_gamma_dtau * P_sb(k,l) * P_sb(m,n) + endif significantShearBandStress + enddo + + endif shearBandingContribution + + call kinetics_twin(Mp,T,dot_gamma_sl,instance,of,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) * f_unrotated + forall (k=1:3,l=1:3,m=1:3,n=1:3) & + dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + + ddot_gamma_dtau_twin(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i) * f_unrotated + enddo twinContibution + + call kinetics_trans(Mp,T,dot_gamma_sl,instance,of,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) * f_unrotated + forall (k=1:3,l=1:3,m=1:3,n=1:3) & + dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + + ddot_gamma_dtau_trans(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i) * f_unrotated + enddo transContibution + + + end associate end subroutine plastic_dislotwin_LpAndItsTangent @@ -707,100 +708,100 @@ end subroutine plastic_dislotwin_LpAndItsTangent !-------------------------------------------------------------------------------------------------- module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) - real(pReal), dimension(3,3), intent(in):: & - Mp !< Mandel stress - real(pReal), intent(in) :: & - T !< temperature at integration point - integer, intent(in) :: & - instance, & - of + real(pReal), dimension(3,3), intent(in):: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + T !< temperature at integration point + integer, intent(in) :: & + instance, & + of - integer :: i - real(pReal) :: & - f_unrotated, & - VacancyDiffusion, & - rho_dip_distance, & - v_cl, & !< climb velocity - Gamma, & !< stacking fault energy - tau, & - sigma_cl, & !< climb stress - b_d !< ratio of burgers vector to stacking fault width - real(pReal), dimension(param(instance)%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) :: & - dot_gamma_twin - real(pReal), dimension(param(instance)%sum_N_tr) :: & - dot_gamma_tr + integer :: i + real(pReal) :: & + f_unrotated, & + VacancyDiffusion, & + rho_dip_distance, & + v_cl, & !< climb velocity + Gamma, & !< stacking fault energy + tau, & + sigma_cl, & !< climb stress + b_d !< ratio of burgers vector to stacking fault width + real(pReal), dimension(param(instance)%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) :: & + dot_gamma_twin + real(pReal), dimension(param(instance)%sum_N_tr) :: & + dot_gamma_tr - associate(prm => param(instance), stt => state(instance), & - dot => dotState(instance), dst => dependentState(instance)) + associate(prm => param(instance), stt => state(instance), & + dot => dotState(instance), dst => dependentState(instance)) - 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)) - VacancyDiffusion = prm%D0*exp(-prm%Qsd/(kB*T)) + 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)) + VacancyDiffusion = prm%D0*exp(-prm%Qsd/(kB*T)) - call kinetics_slip(Mp,T,instance,of,dot_gamma_sl) - dot%gamma_sl(:,of) = abs(dot_gamma_sl) - - rho_dip_distance_min = prm%CEdgeDipMinDistance*prm%b_sl - - slipState: do i = 1, prm%sum_N_sl - tau = math_mul33xx33(Mp,prm%P_sl(1:3,1:3,i)) + call kinetics_slip(Mp,T,instance,of,dot_gamma_sl) + dot%gamma_sl(:,of) = abs(dot_gamma_sl) + + rho_dip_distance_min = prm%CEdgeDipMinDistance*prm%b_sl + + slipState: do i = 1, prm%sum_N_sl + tau = math_mul33xx33(Mp,prm%P_sl(1:3,1:3,i)) - significantSlipStress: if (dEq0(tau)) then - dot_rho_dip_formation(i) = 0.0_pReal - 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, left = rho_dip_distance_min(i)) + significantSlipStress: if (dEq0(tau)) then + dot_rho_dip_formation(i) = 0.0_pReal + 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, 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)) - else - dot_rho_dip_formation(i) = 0.0_pReal - endif + 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)) + else + dot_rho_dip_formation(i) = 0.0_pReal + endif - if (dEq(rho_dip_distance,rho_dip_distance_min(i))) then - dot_rho_dip_climb(i) = 0.0_pReal - else - !@details: Refer: Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981 - sigma_cl = dot_product(prm%n0_sl(1:3,i),matmul(Mp,prm%n0_sl(1:3,i))) - if (prm%ExtendedDislocations) then - Gamma = prm%SFE_0K + prm%dSFE_dT * T - b_d = 24.0_pReal*PI*(1.0_pReal - prm%nu)/(2.0_pReal + prm%nu)* Gamma/(prm%mu*prm%b_sl(i)) - else - b_d = 1.0_pReal - endif - v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Qsd/(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) & - / (rho_dip_distance-rho_dip_distance_min(i)) - endif - endif significantSlipStress - enddo slipState + if (dEq(rho_dip_distance,rho_dip_distance_min(i))) then + dot_rho_dip_climb(i) = 0.0_pReal + else + !@details: Refer: Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981 + sigma_cl = dot_product(prm%n0_sl(1:3,i),matmul(Mp,prm%n0_sl(1:3,i))) + if (prm%ExtendedDislocations) then + Gamma = prm%SFE_0K + prm%dSFE_dT * T + b_d = 24.0_pReal*PI*(1.0_pReal - prm%nu)/(2.0_pReal + prm%nu)* Gamma/(prm%mu*prm%b_sl(i)) + else + b_d = 1.0_pReal + endif + v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Qsd/(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) & + / (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_dip_formation & - - 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_mob(:,of)*abs(dot_gamma_sl) + dot%rho_mob(:,of) = abs(dot_gamma_sl)/(prm%b_sl*dst%Lambda_sl(:,of)) & + - dot_rho_dip_formation & + - 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_mob(:,of)*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_climb + 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_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,instance,of,dot_gamma_twin) + dot%f_tw(:,of) = 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,instance,of,dot_gamma_tr) + dot%f_tr(:,of) = f_unrotated*dot_gamma_tr - end associate + end associate end subroutine plastic_dislotwin_dotState @@ -810,89 +811,89 @@ end subroutine plastic_dislotwin_dotState !-------------------------------------------------------------------------------------------------- module subroutine plastic_dislotwin_dependentState(T,instance,of) - integer, intent(in) :: & - instance, & - of - real(pReal), intent(in) :: & - T - - real(pReal) :: & - sumf_twin,Gamma,sumf_trans - real(pReal), dimension(param(instance)%sum_N_sl) :: & - inv_lambda_sl_sl, & !< 1/mean free distance between 2 forest dislocations seen by a moving dislocation - inv_lambda_sl_tw, & !< 1/mean free distance between 2 twin stacks from different systems seen by a moving dislocation - inv_lambda_sl_tr !< 1/mean free distance between 2 martensite lamellar from different systems seen by a moving dislocation - real(pReal), dimension(param(instance)%sum_N_tw) :: & - inv_lambda_tw_tw, & !< 1/mean free distance between 2 twin stacks from different systems seen by a growing twin - f_over_t_tw - real(pReal), dimension(param(instance)%sum_N_tr) :: & - inv_lambda_tr_tr, & !< 1/mean free distance between 2 martensite stacks from different systems seen by a growing martensite - f_over_t_tr - real(pReal), dimension(:), allocatable :: & - x0 - - - associate(prm => param(instance),& - 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)) - - Gamma = prm%SFE_0K + prm%dSFE_dT * T + integer, intent(in) :: & + instance, & + of + real(pReal), intent(in) :: & + 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_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%CLambdaSlip - - 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) - - inv_lambda_tw_tw = matmul(prm%h_tw_tw,f_over_t_tw)/(1.0_pReal-sumf_twin) + real(pReal) :: & + sumf_twin,Gamma,sumf_trans + real(pReal), dimension(param(instance)%sum_N_sl) :: & + inv_lambda_sl_sl, & !< 1/mean free distance between 2 forest dislocations seen by a moving dislocation + inv_lambda_sl_tw, & !< 1/mean free distance between 2 twin stacks from different systems seen by a moving dislocation + inv_lambda_sl_tr !< 1/mean free distance between 2 martensite lamellar from different systems seen by a moving dislocation + real(pReal), dimension(param(instance)%sum_N_tw) :: & + inv_lambda_tw_tw, & !< 1/mean free distance between 2 twin stacks from different systems seen by a growing twin + f_over_t_tw + real(pReal), dimension(param(instance)%sum_N_tr) :: & + inv_lambda_tr_tr, & !< 1/mean free distance between 2 martensite stacks from different systems seen by a growing martensite + f_over_t_tr + real(pReal), dimension(:), allocatable :: & + x0 - if (prm%sum_N_tr > 0 .and. prm%sum_N_sl > 0) & - inv_lambda_sl_tr = matmul(prm%h_sl_tr,f_over_t_tr)/(1.0_pReal-sumf_trans) - 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 & - / (1.0_pReal+prm%D*(inv_lambda_sl_sl + inv_lambda_sl_tw + inv_lambda_sl_tr)) - else - dst%Lambda_sl(:,of) = 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) - - !* 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))) - - !* 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) & - + 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) & - + 3.0_pReal*prm%b_tr*prm%mu/(prm%L_tr*prm%b_sl) & ! slip burgers here correct? - + prm%h*prm%gamma_fcc_hex/ (3.0_pReal*prm%b_tr) + associate(prm => param(instance),& + stt => state(instance),& + dst => dependentState(instance)) - 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 - - - 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%xc_twin)+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%xc_trans)+cos(pi/3.0_pReal)/x0) - - end associate + sumf_twin = sum(stt%f_tw(1:prm%sum_N_tw,of)) + sumf_trans = sum(stt%f_tr(1:prm%sum_N_tr,of)) + + Gamma = prm%SFE_0K + prm%dSFE_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_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%CLambdaSlip + + 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) + + inv_lambda_tw_tw = matmul(prm%h_tw_tw,f_over_t_tw)/(1.0_pReal-sumf_twin) + + if (prm%sum_N_tr > 0 .and. prm%sum_N_sl > 0) & + inv_lambda_sl_tr = matmul(prm%h_sl_tr,f_over_t_tr)/(1.0_pReal-sumf_trans) + + 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 & + / (1.0_pReal+prm%D*(inv_lambda_sl_sl + inv_lambda_sl_tw + inv_lambda_sl_tr)) + else + dst%Lambda_sl(:,of) = 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) + + !* 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))) + + !* 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) & + + 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) & + + 3.0_pReal*prm%b_tr*prm%mu/(prm%L_tr*prm%b_sl) & ! slip burgers here correct? + + prm%h*prm%gamma_fcc_hex/ (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 + + + 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%xc_twin)+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%xc_trans)+cos(pi/3.0_pReal)/x0) + + end associate end subroutine plastic_dislotwin_dependentState @@ -957,69 +958,69 @@ end subroutine plastic_dislotwin_results pure subroutine kinetics_slip(Mp,T,instance,of, & dot_gamma_sl,ddot_gamma_dtau_slip,tau_slip) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - real(pReal), intent(in) :: & - T !< temperature - integer, intent(in) :: & - instance, & - of - - real(pReal), dimension(param(instance)%sum_N_sl), intent(out) :: & - dot_gamma_sl - real(pReal), dimension(param(instance)%sum_N_sl), optional, intent(out) :: & - ddot_gamma_dtau_slip, & - tau_slip - real(pReal), dimension(param(instance)%sum_N_sl) :: & - ddot_gamma_dtau - - real(pReal), dimension(param(instance)%sum_N_sl) :: & - tau, & - stressRatio, & - StressRatio_p, & - BoltzmannRatio, & - v_wait_inverse, & !< inverse of the effective velocity of a dislocation waiting at obstacles (unsigned) - v_run_inverse, & !< inverse of the velocity of a free moving dislocation (unsigned) - dV_wait_inverse_dTau, & - dV_run_inverse_dTau, & - dV_dTau, & - tau_eff !< effective resolved stress - integer :: i + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + T !< temperature + integer, intent(in) :: & + instance, & + of + + real(pReal), dimension(param(instance)%sum_N_sl), intent(out) :: & + dot_gamma_sl + real(pReal), dimension(param(instance)%sum_N_sl), optional, intent(out) :: & + ddot_gamma_dtau_slip, & + tau_slip + real(pReal), dimension(param(instance)%sum_N_sl) :: & + ddot_gamma_dtau - associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) - - do i = 1, prm%sum_N_sl - tau(i) = math_mul33xx33(Mp,prm%P_sl(1:3,1:3,i)) - enddo - - tau_eff = abs(tau)-dst%tau_pass(:,of) - - significantStress: where(tau_eff > tol_math_check) - stressRatio = tau_eff/prm%tau_0 - StressRatio_p = stressRatio** prm%p - BoltzmannRatio = prm%Delta_F/(kB*T) - v_wait_inverse = prm%v0**(-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) - - dV_wait_inverse_dTau = -1.0_pReal * v_wait_inverse * prm%p * prm%q * BoltzmannRatio & - * (stressRatio**(prm%p-1.0_pReal)) & - * (1.0_pReal-StressRatio_p)**(prm%q-1.0_pReal) & - / prm%tau_0 - 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 - else where significantStress - dot_gamma_sl = 0.0_pReal - ddot_gamma_dtau = 0.0_pReal - end where significantStress - - end associate + real(pReal), dimension(param(instance)%sum_N_sl) :: & + tau, & + stressRatio, & + StressRatio_p, & + BoltzmannRatio, & + v_wait_inverse, & !< inverse of the effective velocity of a dislocation waiting at obstacles (unsigned) + v_run_inverse, & !< inverse of the velocity of a free moving dislocation (unsigned) + dV_wait_inverse_dTau, & + dV_run_inverse_dTau, & + dV_dTau, & + tau_eff !< effective resolved stress + integer :: i - if(present(ddot_gamma_dtau_slip)) ddot_gamma_dtau_slip = ddot_gamma_dtau - if(present(tau_slip)) tau_slip = tau + associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) + + do i = 1, prm%sum_N_sl + tau(i) = math_mul33xx33(Mp,prm%P_sl(1:3,1:3,i)) + enddo + + tau_eff = abs(tau)-dst%tau_pass(:,of) + + significantStress: where(tau_eff > tol_math_check) + stressRatio = tau_eff/prm%tau_0 + StressRatio_p = stressRatio** prm%p + BoltzmannRatio = prm%Delta_F/(kB*T) + v_wait_inverse = prm%v0**(-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) + + dV_wait_inverse_dTau = -1.0_pReal * v_wait_inverse * prm%p * prm%q * BoltzmannRatio & + * (stressRatio**(prm%p-1.0_pReal)) & + * (1.0_pReal-StressRatio_p)**(prm%q-1.0_pReal) & + / prm%tau_0 + 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 + else where significantStress + dot_gamma_sl = 0.0_pReal + ddot_gamma_dtau = 0.0_pReal + end where significantStress + + end associate + + if(present(ddot_gamma_dtau_slip)) ddot_gamma_dtau_slip = ddot_gamma_dtau + if(present(tau_slip)) tau_slip = tau end subroutine kinetics_slip @@ -1030,61 +1031,61 @@ end subroutine kinetics_slip pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,of,& dot_gamma_twin,ddot_gamma_dtau_twin) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - real(pReal), intent(in) :: & - T !< temperature - integer, intent(in) :: & - instance, & - of - real(pReal), dimension(param(instance)%sum_N_sl), intent(in) :: & - dot_gamma_sl - - real(pReal), dimension(param(instance)%sum_N_tw), intent(out) :: & - dot_gamma_twin - real(pReal), dimension(param(instance)%sum_N_tw), optional, intent(out) :: & - ddot_gamma_dtau_twin - - real, dimension(param(instance)%sum_N_tw) :: & - tau, & - Ndot0, & - stressRatio_r, & - ddot_gamma_dtau - - integer :: i,s1,s2 + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + T !< temperature + integer, intent(in) :: & + instance, & + of + real(pReal), dimension(param(instance)%sum_N_sl), intent(in) :: & + dot_gamma_sl + + real(pReal), dimension(param(instance)%sum_N_tw), intent(out) :: & + dot_gamma_twin + real(pReal), dimension(param(instance)%sum_N_tw), optional, intent(out) :: & + ddot_gamma_dtau_twin - associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) - - do i = 1, prm%sum_N_tw - tau(i) = math_mul33xx33(Mp,prm%P_tw(1:3,1:3,i)) - 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 - (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 - else - Ndot0=0.0_pReal - end if - else isFCC - Ndot0=prm%dot_N_0_tw(i) - endif isFCC - 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) - ddot_gamma_dtau = (dot_gamma_twin*prm%r/tau)*StressRatio_r - else where significantStress - dot_gamma_twin = 0.0_pReal - ddot_gamma_dtau = 0.0_pReal - end where significantStress + real, dimension(param(instance)%sum_N_tw) :: & + tau, & + Ndot0, & + stressRatio_r, & + ddot_gamma_dtau - end associate - - if(present(ddot_gamma_dtau_twin)) ddot_gamma_dtau_twin = ddot_gamma_dtau + integer :: i,s1,s2 + + associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) + + do i = 1, prm%sum_N_tw + tau(i) = math_mul33xx33(Mp,prm%P_tw(1:3,1:3,i)) + 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 + (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 + else + Ndot0=0.0_pReal + end if + else isFCC + Ndot0=prm%dot_N_0_tw(i) + endif isFCC + 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) + ddot_gamma_dtau = (dot_gamma_twin*prm%r/tau)*StressRatio_r + else where significantStress + dot_gamma_twin = 0.0_pReal + ddot_gamma_dtau = 0.0_pReal + end where significantStress + + end associate + + if(present(ddot_gamma_dtau_twin)) ddot_gamma_dtau_twin = ddot_gamma_dtau end subroutine kinetics_twin @@ -1095,60 +1096,60 @@ end subroutine kinetics_twin pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,of,& dot_gamma_tr,ddot_gamma_dtau_trans) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - real(pReal), intent(in) :: & - T !< temperature - integer, intent(in) :: & - instance, & - of - real(pReal), dimension(param(instance)%sum_N_sl), intent(in) :: & - dot_gamma_sl - - real(pReal), dimension(param(instance)%sum_N_tr), intent(out) :: & - dot_gamma_tr - real(pReal), dimension(param(instance)%sum_N_tr), optional, intent(out) :: & - ddot_gamma_dtau_trans - - real, dimension(param(instance)%sum_N_tr) :: & - tau, & - Ndot0, & - stressRatio_s, & - ddot_gamma_dtau - - integer :: i,s1,s2 - associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) - - do i = 1, prm%sum_N_tr - tau(i) = math_mul33xx33(Mp,prm%P_tr(1:3,1:3,i)) - 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 - (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 - else - Ndot0=0.0_pReal - end if - else isFCC - Ndot0=prm%dot_N_0_tr(i) - endif isFCC - 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) - ddot_gamma_dtau = (dot_gamma_tr*prm%s/tau)*StressRatio_s - else where significantStress - dot_gamma_tr = 0.0_pReal - ddot_gamma_dtau = 0.0_pReal - end where significantStress + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + T !< temperature + integer, intent(in) :: & + instance, & + of + real(pReal), dimension(param(instance)%sum_N_sl), intent(in) :: & + dot_gamma_sl + + real(pReal), dimension(param(instance)%sum_N_tr), intent(out) :: & + dot_gamma_tr + real(pReal), dimension(param(instance)%sum_N_tr), optional, intent(out) :: & + ddot_gamma_dtau_trans - end associate - - if(present(ddot_gamma_dtau_trans)) ddot_gamma_dtau_trans = ddot_gamma_dtau + real, dimension(param(instance)%sum_N_tr) :: & + tau, & + Ndot0, & + stressRatio_s, & + ddot_gamma_dtau + + integer :: i,s1,s2 + associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) + + do i = 1, prm%sum_N_tr + tau(i) = math_mul33xx33(Mp,prm%P_tr(1:3,1:3,i)) + 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 + (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 + else + Ndot0=0.0_pReal + end if + else isFCC + Ndot0=prm%dot_N_0_tr(i) + endif isFCC + 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) + ddot_gamma_dtau = (dot_gamma_tr*prm%s/tau)*StressRatio_s + else where significantStress + dot_gamma_tr = 0.0_pReal + ddot_gamma_dtau = 0.0_pReal + end where significantStress + + end associate + + if(present(ddot_gamma_dtau_trans)) ddot_gamma_dtau_trans = ddot_gamma_dtau end subroutine kinetics_trans diff --git a/src/constitutive_plastic_kinehardening.f90 b/src/constitutive_plastic_kinehardening.f90 index 2957ff80c..4949c6533 100644 --- a/src/constitutive_plastic_kinehardening.f90 +++ b/src/constitutive_plastic_kinehardening.f90 @@ -7,66 +7,66 @@ !-------------------------------------------------------------------------------------------------- submodule(constitutive) plastic_kinehardening - enum, bind(c) - enumerator :: & - undefined_ID, & - crss_ID, & !< critical resolved stress - crss_back_ID, & !< critical resolved back stress - sense_ID, & !< sense of acting shear stress (-1 or +1) - chi0_ID, & !< backstress at last switch of stress sense (positive?) - gamma0_ID, & !< accumulated shear at last switch of stress sense (at current switch?) - accshear_ID, & - shearrate_ID, & - resolvedstress_ID - end enum - - type :: tParameters - real(pReal) :: & - gdot0, & !< reference shear strain rate for slip - n, & !< stress exponent for slip - aTolResistance, & - aTolShear - real(pReal), allocatable, dimension(:) :: & - crss0, & !< initial critical shear stress for slip - theta0, & !< initial hardening rate of forward stress for each slip - theta1, & !< asymptotic hardening rate of forward stress for each slip - theta0_b, & !< initial hardening rate of back stress for each slip - theta1_b, & !< asymptotic hardening rate of back stress for each slip - tau1, & - tau1_b, & - nonSchmidCoeff - real(pReal), allocatable, dimension(:,:) :: & - interaction_slipslip !< slip resistance from slip activity - real(pReal), allocatable, dimension(:,:,:) :: & - Schmid, & - nonSchmid_pos, & - nonSchmid_neg - integer :: & - totalNslip, & !< total number of active slip system - of_debug = 0 - integer, allocatable, dimension(:) :: & - Nslip !< number of active slip systems for each family - integer(kind(undefined_ID)), allocatable, dimension(:) :: & - outputID !< ID of each post result output - end type tParameters - - type :: tKinehardeningState - real(pReal), pointer, dimension(:,:) :: & !< vectors along NipcMyInstance - crss, & !< critical resolved stress - crss_back, & !< critical resolved back stress - sense, & !< sense of acting shear stress (-1 or +1) - chi0, & !< backstress at last switch of stress sense - gamma0, & !< accumulated shear at last switch of stress sense - accshear !< accumulated (absolute) shear - end type tKinehardeningState - + enum, bind(c) + enumerator :: & + undefined_ID, & + crss_ID, & !< critical resolved stress + crss_back_ID, & !< critical resolved back stress + sense_ID, & !< sense of acting shear stress (-1 or +1) + chi0_ID, & !< backstress at last switch of stress sense (positive?) + gamma0_ID, & !< accumulated shear at last switch of stress sense (at current switch?) + accshear_ID, & + shearrate_ID, & + resolvedstress_ID + end enum + + type :: tParameters + real(pReal) :: & + gdot0, & !< reference shear strain rate for slip + n, & !< stress exponent for slip + aTolResistance, & + aTolShear + real(pReal), allocatable, dimension(:) :: & + crss0, & !< initial critical shear stress for slip + theta0, & !< initial hardening rate of forward stress for each slip + theta1, & !< asymptotic hardening rate of forward stress for each slip + theta0_b, & !< initial hardening rate of back stress for each slip + theta1_b, & !< asymptotic hardening rate of back stress for each slip + tau1, & + tau1_b, & + nonSchmidCoeff + real(pReal), allocatable, dimension(:,:) :: & + interaction_slipslip !< slip resistance from slip activity + real(pReal), allocatable, dimension(:,:,:) :: & + Schmid, & + nonSchmid_pos, & + nonSchmid_neg + integer :: & + totalNslip, & !< total number of active slip system + of_debug = 0 + integer, allocatable, dimension(:) :: & + Nslip !< number of active slip systems for each family + integer(kind(undefined_ID)), allocatable, dimension(:) :: & + outputID !< ID of each post result output + end type tParameters + + type :: tKinehardeningState + real(pReal), pointer, dimension(:,:) :: & !< vectors along NipcMyInstance + crss, & !< critical resolved stress + crss_back, & !< critical resolved back stress + sense, & !< sense of acting shear stress (-1 or +1) + chi0, & !< backstress at last switch of stress sense + gamma0, & !< accumulated shear at last switch of stress sense + accshear !< accumulated (absolute) shear + end type tKinehardeningState + !-------------------------------------------------------------------------------------------------- ! containers for parameters and state - type(tParameters), allocatable, dimension(:) :: param - type(tKinehardeningState), allocatable, dimension(:) :: & - dotState, & - deltaState, & - state + type(tParameters), allocatable, dimension(:) :: param + type(tKinehardeningState), allocatable, dimension(:) :: & + dotState, & + deltaState, & + state contains @@ -77,202 +77,202 @@ contains !-------------------------------------------------------------------------------------------------- module subroutine plastic_kinehardening_init - integer :: & - Ninstance, & - p, i, o, & - NipcMyPhase, & - sizeState, sizeDeltaState, sizeDotState, & - startIndex, endIndex - - integer(kind(undefined_ID)) :: & - outputID - - character(len=pStringLen) :: & - extmsg = '' - character(len=pStringLen), dimension(:), allocatable :: & - outputs - - write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_KINEHARDENING_label//' init -+>>>' - - Ninstance = count(phase_plasticity == PLASTICITY_KINEHARDENING_ID) - if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & - write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - - allocate(param(Ninstance)) - allocate(state(Ninstance)) - allocate(dotState(Ninstance)) - allocate(deltaState(Ninstance)) - - do p = 1, size(phase_plasticityInstance) - if (phase_plasticity(p) /= PLASTICITY_KINEHARDENING_ID) cycle - associate(prm => param(phase_plasticityInstance(p)), & - dot => dotState(phase_plasticityInstance(p)), & - dlt => deltaState(phase_plasticityInstance(p)), & - stt => state(phase_plasticityInstance(p)),& - config => config_phase(p)) - + integer :: & + Ninstance, & + p, i, o, & + NipcMyPhase, & + sizeState, sizeDeltaState, sizeDotState, & + startIndex, endIndex + + integer(kind(undefined_ID)) :: & + outputID + + character(len=pStringLen) :: & + extmsg = '' + character(len=pStringLen), dimension(:), allocatable :: & + outputs + + write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_KINEHARDENING_label//' init -+>>>' + + Ninstance = count(phase_plasticity == PLASTICITY_KINEHARDENING_ID) + if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & + write(6,'(a16,1x,i5,/)') '# instances:',Ninstance + + allocate(param(Ninstance)) + allocate(state(Ninstance)) + allocate(dotState(Ninstance)) + allocate(deltaState(Ninstance)) + + do p = 1, size(phase_plasticityInstance) + if (phase_plasticity(p) /= PLASTICITY_KINEHARDENING_ID) cycle + associate(prm => param(phase_plasticityInstance(p)), & + dot => dotState(phase_plasticityInstance(p)), & + dlt => deltaState(phase_plasticityInstance(p)), & + stt => state(phase_plasticityInstance(p)),& + config => config_phase(p)) + #ifdef DEBUG - if (p==material_phaseAt(debug_g,debug_e)) then - prm%of_debug = material_phasememberAt(debug_g,debug_i,debug_e) - endif + if (p==material_phaseAt(debug_g,debug_e)) then + prm%of_debug = material_phasememberAt(debug_g,debug_i,debug_e) + endif #endif - + !-------------------------------------------------------------------------------------------------- ! optional parameters that need to be defined - prm%aTolResistance = config%getFloat('atol_resistance',defaultVal=1.0_pReal) - prm%aTolShear = config%getFloat('atol_shear', defaultVal=1.0e-6_pReal) - - ! sanity checks - if (prm%aTolResistance <= 0.0_pReal) extmsg = trim(extmsg)//' aTolresistance' - if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//' aTolShear' - + prm%aTolResistance = config%getFloat('atol_resistance',defaultVal=1.0_pReal) + prm%aTolShear = config%getFloat('atol_shear', defaultVal=1.0e-6_pReal) + + ! sanity checks + if (prm%aTolResistance <= 0.0_pReal) extmsg = trim(extmsg)//' aTolresistance' + if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//' aTolShear' + !-------------------------------------------------------------------------------------------------- ! slip related parameters - prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) - prm%totalNslip = sum(prm%Nslip) - slipActive: if (prm%totalNslip > 0) then - prm%Schmid = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),& - config%getFloat('c/a',defaultVal=0.0_pReal)) - - if(trim(config%getString('lattice_structure')) == 'bcc') then - prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',& - defaultVal = emptyRealArray) - prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1) - prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1) - else - prm%nonSchmid_pos = prm%Schmid - prm%nonSchmid_neg = prm%Schmid - endif - prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(prm%Nslip, & - config%getFloats('interaction_slipslip'), & - config%getString('lattice_structure')) - - prm%crss0 = config%getFloats('crss0', requiredSize=size(prm%Nslip)) - prm%tau1 = config%getFloats('tau1', requiredSize=size(prm%Nslip)) - prm%tau1_b = config%getFloats('tau1_b', requiredSize=size(prm%Nslip)) - prm%theta0 = config%getFloats('theta0', requiredSize=size(prm%Nslip)) - prm%theta1 = config%getFloats('theta1', requiredSize=size(prm%Nslip)) - prm%theta0_b = config%getFloats('theta0_b', requiredSize=size(prm%Nslip)) - prm%theta1_b = config%getFloats('theta1_b', requiredSize=size(prm%Nslip)) - - prm%gdot0 = config%getFloat('gdot0') - prm%n = config%getFloat('n_slip') - - ! expand: family => system - prm%crss0 = math_expand(prm%crss0, prm%Nslip) - prm%tau1 = math_expand(prm%tau1, prm%Nslip) - prm%tau1_b = math_expand(prm%tau1_b, prm%Nslip) - prm%theta0 = math_expand(prm%theta0, prm%Nslip) - prm%theta1 = math_expand(prm%theta1, prm%Nslip) - prm%theta0_b = math_expand(prm%theta0_b,prm%Nslip) - prm%theta1_b = math_expand(prm%theta1_b,prm%Nslip) - - - + prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) + prm%totalNslip = sum(prm%Nslip) + slipActive: if (prm%totalNslip > 0) then + prm%Schmid = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),& + config%getFloat('c/a',defaultVal=0.0_pReal)) + + if(trim(config%getString('lattice_structure')) == 'bcc') then + prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',& + defaultVal = emptyRealArray) + prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1) + prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1) + else + prm%nonSchmid_pos = prm%Schmid + prm%nonSchmid_neg = prm%Schmid + endif + prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(prm%Nslip, & + config%getFloats('interaction_slipslip'), & + config%getString('lattice_structure')) + + prm%crss0 = config%getFloats('crss0', requiredSize=size(prm%Nslip)) + prm%tau1 = config%getFloats('tau1', requiredSize=size(prm%Nslip)) + prm%tau1_b = config%getFloats('tau1_b', requiredSize=size(prm%Nslip)) + prm%theta0 = config%getFloats('theta0', requiredSize=size(prm%Nslip)) + prm%theta1 = config%getFloats('theta1', requiredSize=size(prm%Nslip)) + prm%theta0_b = config%getFloats('theta0_b', requiredSize=size(prm%Nslip)) + prm%theta1_b = config%getFloats('theta1_b', requiredSize=size(prm%Nslip)) + + prm%gdot0 = config%getFloat('gdot0') + prm%n = config%getFloat('n_slip') + + ! expand: family => system + prm%crss0 = math_expand(prm%crss0, prm%Nslip) + prm%tau1 = math_expand(prm%tau1, prm%Nslip) + prm%tau1_b = math_expand(prm%tau1_b, prm%Nslip) + prm%theta0 = math_expand(prm%theta0, prm%Nslip) + prm%theta1 = math_expand(prm%theta1, prm%Nslip) + prm%theta0_b = math_expand(prm%theta0_b,prm%Nslip) + prm%theta1_b = math_expand(prm%theta1_b,prm%Nslip) + + + !-------------------------------------------------------------------------------------------------- ! sanity checks - if ( prm%gdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0' - if ( prm%n <= 0.0_pReal) extmsg = trim(extmsg)//' n_slip' - if (any(prm%crss0 <= 0.0_pReal)) extmsg = trim(extmsg)//' crss0' - if (any(prm%tau1 <= 0.0_pReal)) extmsg = trim(extmsg)//' tau1' - if (any(prm%tau1_b <= 0.0_pReal)) extmsg = trim(extmsg)//' tau1_b' - - !ToDo: Any sensible checks for theta? - - endif slipActive - + if ( prm%gdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0' + if ( prm%n <= 0.0_pReal) extmsg = trim(extmsg)//' n_slip' + if (any(prm%crss0 <= 0.0_pReal)) extmsg = trim(extmsg)//' crss0' + if (any(prm%tau1 <= 0.0_pReal)) extmsg = trim(extmsg)//' tau1' + if (any(prm%tau1_b <= 0.0_pReal)) extmsg = trim(extmsg)//' tau1_b' + + !ToDo: Any sensible checks for theta? + + endif slipActive + !-------------------------------------------------------------------------------------------------- ! exit if any parameter is out of range - if (extmsg /= '') & - call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_KINEHARDENING_label//')') - + if (extmsg /= '') & + call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_KINEHARDENING_label//')') + !-------------------------------------------------------------------------------------------------- ! output pararameters - outputs = config%getStrings('(output)',defaultVal=emptyStringArray) - allocate(prm%outputID(0)) - do i=1, size(outputs) - outputID = undefined_ID - select case(outputs(i)) - - case ('resistance') - outputID = merge(crss_ID,undefined_ID,prm%totalNslip>0) - case ('accumulatedshear') - outputID = merge(accshear_ID,undefined_ID,prm%totalNslip>0) - case ('shearrate') - outputID = merge(shearrate_ID,undefined_ID,prm%totalNslip>0) - case ('resolvedstress') - outputID = merge(resolvedstress_ID,undefined_ID,prm%totalNslip>0) - case ('backstress') - outputID = merge(crss_back_ID,undefined_ID,prm%totalNslip>0) - case ('sense') - outputID = merge(sense_ID,undefined_ID,prm%totalNslip>0) - case ('chi0') - outputID = merge(chi0_ID,undefined_ID,prm%totalNslip>0) - case ('gamma0') - outputID = merge(gamma0_ID,undefined_ID,prm%totalNslip>0) - - end select - - if (outputID /= undefined_ID) then - prm%outputID = [prm%outputID , outputID] - endif - - enddo - + outputs = config%getStrings('(output)',defaultVal=emptyStringArray) + allocate(prm%outputID(0)) + do i=1, size(outputs) + outputID = undefined_ID + select case(outputs(i)) + + case ('resistance') + outputID = merge(crss_ID,undefined_ID,prm%totalNslip>0) + case ('accumulatedshear') + outputID = merge(accshear_ID,undefined_ID,prm%totalNslip>0) + case ('shearrate') + outputID = merge(shearrate_ID,undefined_ID,prm%totalNslip>0) + case ('resolvedstress') + outputID = merge(resolvedstress_ID,undefined_ID,prm%totalNslip>0) + case ('backstress') + outputID = merge(crss_back_ID,undefined_ID,prm%totalNslip>0) + case ('sense') + outputID = merge(sense_ID,undefined_ID,prm%totalNslip>0) + case ('chi0') + outputID = merge(chi0_ID,undefined_ID,prm%totalNslip>0) + case ('gamma0') + outputID = merge(gamma0_ID,undefined_ID,prm%totalNslip>0) + + end select + + if (outputID /= undefined_ID) then + prm%outputID = [prm%outputID , outputID] + endif + + enddo + !-------------------------------------------------------------------------------------------------- ! allocate state arrays - NipcMyPhase = count(material_phaseAt == p) * discretization_nIP - sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%totalNslip - sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%totalNslip - sizeState = sizeDotState + sizeDeltaState - - call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,sizeDeltaState) - + NipcMyPhase = count(material_phaseAt == p) * discretization_nIP + sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%totalNslip + sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%totalNslip + sizeState = sizeDotState + sizeDeltaState + + call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,sizeDeltaState) + !-------------------------------------------------------------------------------------------------- ! locally defined state aliases and initialization of state0 and aTolState - startIndex = 1 - endIndex = prm%totalNslip - stt%crss => plasticState(p)%state (startIndex:endIndex,:) - stt%crss = spread(prm%crss0, 2, NipcMyPhase) - dot%crss => plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance - - startIndex = endIndex + 1 - endIndex = endIndex + prm%totalNslip - stt%crss_back => plasticState(p)%state (startIndex:endIndex,:) - dot%crss_back => plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance - - startIndex = endIndex + 1 - endIndex = endIndex + prm%totalNslip - stt%accshear => plasticState(p)%state (startIndex:endIndex,:) - dot%accshear => plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear - ! global alias - plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%accumulatedSlip => plasticState(p)%state(startIndex:endIndex,:) - - o = plasticState(p)%offsetDeltaState - startIndex = endIndex + 1 - endIndex = endIndex + prm%totalNslip - stt%sense => plasticState(p)%state (startIndex :endIndex ,:) - dlt%sense => plasticState(p)%deltaState(startIndex-o:endIndex-o,:) - - startIndex = endIndex + 1 - endIndex = endIndex + prm%totalNslip - stt%chi0 => plasticState(p)%state (startIndex :endIndex ,:) - dlt%chi0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:) - - startIndex = endIndex + 1 - endIndex = endIndex + prm%totalNslip - stt%gamma0 => plasticState(p)%state (startIndex :endIndex ,:) - dlt%gamma0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:) - - plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally - - end associate - - enddo + startIndex = 1 + endIndex = prm%totalNslip + stt%crss => plasticState(p)%state (startIndex:endIndex,:) + stt%crss = spread(prm%crss0, 2, NipcMyPhase) + dot%crss => plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance + + startIndex = endIndex + 1 + endIndex = endIndex + prm%totalNslip + stt%crss_back => plasticState(p)%state (startIndex:endIndex,:) + dot%crss_back => plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance + + startIndex = endIndex + 1 + endIndex = endIndex + prm%totalNslip + stt%accshear => plasticState(p)%state (startIndex:endIndex,:) + dot%accshear => plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear + ! global alias + plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%accumulatedSlip => plasticState(p)%state(startIndex:endIndex,:) + + o = plasticState(p)%offsetDeltaState + startIndex = endIndex + 1 + endIndex = endIndex + prm%totalNslip + stt%sense => plasticState(p)%state (startIndex :endIndex ,:) + dlt%sense => plasticState(p)%deltaState(startIndex-o:endIndex-o,:) + + startIndex = endIndex + 1 + endIndex = endIndex + prm%totalNslip + stt%chi0 => plasticState(p)%state (startIndex :endIndex ,:) + dlt%chi0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:) + + startIndex = endIndex + 1 + endIndex = endIndex + prm%totalNslip + stt%gamma0 => plasticState(p)%state (startIndex :endIndex ,:) + dlt%gamma0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:) + + plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally + + end associate + + enddo end subroutine plastic_kinehardening_init @@ -282,39 +282,39 @@ end subroutine plastic_kinehardening_init !-------------------------------------------------------------------------------------------------- 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 - - integer :: & - i,k,l,m,n - real(pReal), dimension(param(instance)%totalNslip) :: & - gdot_pos,gdot_neg, & - dgdot_dtau_pos,dgdot_dtau_neg - - Lp = 0.0_pReal - dLp_dMp = 0.0_pReal - - associate(prm => param(instance)) - - call kinetics(Mp,instance,of,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) - - do i = 1, prm%totalNslip - Lp = Lp + (gdot_pos(i)+gdot_neg(i))*prm%Schmid(1:3,1:3,i) - forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & - + dgdot_dtau_pos(i) * prm%Schmid(k,l,i) * prm%nonSchmid_pos(m,n,i) & - + dgdot_dtau_neg(i) * prm%Schmid(k,l,i) * prm%nonSchmid_neg(m,n,i) - enddo - - end associate + 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 + + integer :: & + i,k,l,m,n + real(pReal), dimension(param(instance)%totalNslip) :: & + gdot_pos,gdot_neg, & + dgdot_dtau_pos,dgdot_dtau_neg + + Lp = 0.0_pReal + dLp_dMp = 0.0_pReal + + associate(prm => param(instance)) + + call kinetics(Mp,instance,of,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) + + do i = 1, prm%totalNslip + Lp = Lp + (gdot_pos(i)+gdot_neg(i))*prm%Schmid(1:3,1:3,i) + forall (k=1:3,l=1:3,m=1:3,n=1:3) & + dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + + dgdot_dtau_pos(i) * prm%Schmid(k,l,i) * prm%nonSchmid_pos(m,n,i) & + + dgdot_dtau_neg(i) * prm%Schmid(k,l,i) * prm%nonSchmid_neg(m,n,i) + enddo + + end associate end subroutine plastic_kinehardening_LpAndItsTangent @@ -324,39 +324,39 @@ end subroutine plastic_kinehardening_LpAndItsTangent !-------------------------------------------------------------------------------------------------- module subroutine plastic_kinehardening_dotState(Mp,instance,of) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - - real(pReal) :: & - sumGamma - real(pReal), dimension(param(instance)%totalNslip) :: & - gdot_pos,gdot_neg - - - associate(prm => param(instance), stt => state(instance), dot => dotState(instance)) - - call kinetics(Mp,instance,of,gdot_pos,gdot_neg) - dot%accshear(:,of) = abs(gdot_pos+gdot_neg) - sumGamma = sum(stt%accshear(:,of)) - + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of - dot%crss(:,of) = matmul(prm%interaction_SlipSlip,dot%accshear(:,of)) & - * ( prm%theta1 & - + (prm%theta0 - prm%theta1 + prm%theta0*prm%theta1*sumGamma/prm%tau1) & - * exp(-sumGamma*prm%theta0/prm%tau1) & - ) - - dot%crss_back(:,of) = stt%sense(:,of)*dot%accshear(:,of) * & - ( prm%theta1_b + & - (prm%theta0_b - prm%theta1_b & - + prm%theta0_b*prm%theta1_b/(prm%tau1_b+stt%chi0(:,of))*(stt%accshear(:,of)-stt%gamma0(:,of))& - ) *exp(-(stt%accshear(:,of)-stt%gamma0(:,of)) *prm%theta0_b/(prm%tau1_b+stt%chi0(:,of))) & - ) - - end associate + real(pReal) :: & + sumGamma + real(pReal), dimension(param(instance)%totalNslip) :: & + gdot_pos,gdot_neg + + + associate(prm => param(instance), stt => state(instance), dot => dotState(instance)) + + call kinetics(Mp,instance,of,gdot_pos,gdot_neg) + dot%accshear(:,of) = abs(gdot_pos+gdot_neg) + sumGamma = sum(stt%accshear(:,of)) + + + dot%crss(:,of) = matmul(prm%interaction_SlipSlip,dot%accshear(:,of)) & + * ( prm%theta1 & + + (prm%theta0 - prm%theta1 + prm%theta0*prm%theta1*sumGamma/prm%tau1) & + * exp(-sumGamma*prm%theta0/prm%tau1) & + ) + + dot%crss_back(:,of) = stt%sense(:,of)*dot%accshear(:,of) * & + ( prm%theta1_b + & + (prm%theta0_b - prm%theta1_b & + + prm%theta0_b*prm%theta1_b/(prm%tau1_b+stt%chi0(:,of))*(stt%accshear(:,of)-stt%gamma0(:,of))& + ) *exp(-(stt%accshear(:,of)-stt%gamma0(:,of)) *prm%theta0_b/(prm%tau1_b+stt%chi0(:,of))) & + ) + + end associate end subroutine plastic_kinehardening_dotState @@ -366,45 +366,45 @@ end subroutine plastic_kinehardening_dotState !-------------------------------------------------------------------------------------------------- module subroutine plastic_kinehardening_deltaState(Mp,instance,of) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - - real(pReal), dimension(param(instance)%totalNslip) :: & - gdot_pos,gdot_neg, & - sense - - 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... - 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 - + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + + real(pReal), dimension(param(instance)%totalNslip) :: & + gdot_pos,gdot_neg, & + sense + + 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... + 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 (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0 & - .and. (of == prm%of_debug & - .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0)) then - write(6,'(a)') '======= kinehardening delta state =======' - write(6,*) sense,state(instance)%sense(:,of) - endif + if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0 & + .and. (of == prm%of_debug & + .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0)) then + write(6,'(a)') '======= kinehardening delta state =======' + write(6,*) sense,state(instance)%sense(:,of) + 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 - else where - dlt%sense (:,of) = 0.0_pReal - dlt%chi0 (:,of) = 0.0_pReal - dlt%gamma0(:,of) = 0.0_pReal - end where - - end associate + 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 + else where + dlt%sense (:,of) = 0.0_pReal + dlt%chi0 (:,of) = 0.0_pReal + dlt%gamma0(:,of) = 0.0_pReal + end where + + end associate end subroutine plastic_kinehardening_deltaState @@ -458,64 +458,64 @@ end subroutine plastic_kinehardening_results pure subroutine kinetics(Mp,instance,of, & 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 - - real(pReal), intent(out), dimension(param(instance)%totalNslip) :: & - gdot_pos, & - gdot_neg - real(pReal), intent(out), optional, dimension(param(instance)%totalNslip) :: & - dgdot_dtau_pos, & - dgdot_dtau_neg - - real(pReal), dimension(param(instance)%totalNslip) :: & - tau_pos, & - tau_neg - integer :: i - logical :: nonSchmidActive - - associate(prm => param(instance), stt => state(instance)) - - nonSchmidActive = size(prm%nonSchmidCoeff) > 0 - - do i = 1, prm%totalNslip - tau_pos(i) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i)) - stt%crss_back(i,of) - tau_neg(i) = merge(math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i)) - stt%crss_back(i,of), & - 0.0_pReal, nonSchmidActive) - enddo - - where(dNeq0(tau_pos)) - gdot_pos = prm%gdot0 * merge(0.5_pReal,1.0_pReal, nonSchmidActive) & ! 1/2 if non-Schmid active - * sign(abs(tau_pos/stt%crss(:,of))**prm%n, tau_pos) - else where - gdot_pos = 0.0_pReal - end where - - where(dNeq0(tau_neg)) - gdot_neg = prm%gdot0 * 0.5_pReal & ! only used if non-Schmid active, always 1/2 - * sign(abs(tau_neg/stt%crss(:,of))**prm%n, tau_neg) - else where - gdot_neg = 0.0_pReal - end where - - if (present(dgdot_dtau_pos)) then - where(dNeq0(gdot_pos)) - dgdot_dtau_pos = gdot_pos*prm%n/tau_pos - else where - dgdot_dtau_pos = 0.0_pReal - end where - endif - if (present(dgdot_dtau_neg)) then - where(dNeq0(gdot_neg)) - dgdot_dtau_neg = gdot_neg*prm%n/tau_neg - else where - dgdot_dtau_neg = 0.0_pReal - end where - endif - end associate + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + + real(pReal), intent(out), dimension(param(instance)%totalNslip) :: & + gdot_pos, & + gdot_neg + real(pReal), intent(out), optional, dimension(param(instance)%totalNslip) :: & + dgdot_dtau_pos, & + dgdot_dtau_neg + + real(pReal), dimension(param(instance)%totalNslip) :: & + tau_pos, & + tau_neg + integer :: i + logical :: nonSchmidActive + + associate(prm => param(instance), stt => state(instance)) + + nonSchmidActive = size(prm%nonSchmidCoeff) > 0 + + do i = 1, prm%totalNslip + tau_pos(i) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i)) - stt%crss_back(i,of) + tau_neg(i) = merge(math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i)) - stt%crss_back(i,of), & + 0.0_pReal, nonSchmidActive) + enddo + + where(dNeq0(tau_pos)) + gdot_pos = prm%gdot0 * merge(0.5_pReal,1.0_pReal, nonSchmidActive) & ! 1/2 if non-Schmid active + * sign(abs(tau_pos/stt%crss(:,of))**prm%n, tau_pos) + else where + gdot_pos = 0.0_pReal + end where + + where(dNeq0(tau_neg)) + gdot_neg = prm%gdot0 * 0.5_pReal & ! only used if non-Schmid active, always 1/2 + * sign(abs(tau_neg/stt%crss(:,of))**prm%n, tau_neg) + else where + gdot_neg = 0.0_pReal + end where + + if (present(dgdot_dtau_pos)) then + where(dNeq0(gdot_pos)) + dgdot_dtau_pos = gdot_pos*prm%n/tau_pos + else where + dgdot_dtau_pos = 0.0_pReal + end where + endif + if (present(dgdot_dtau_neg)) then + where(dNeq0(gdot_neg)) + dgdot_dtau_neg = gdot_neg*prm%n/tau_neg + else where + dgdot_dtau_neg = 0.0_pReal + end where + endif + end associate end subroutine kinetics diff --git a/src/constitutive_plastic_none.f90 b/src/constitutive_plastic_none.f90 index 5c336e2d7..578fb4149 100644 --- a/src/constitutive_plastic_none.f90 +++ b/src/constitutive_plastic_none.f90 @@ -14,24 +14,24 @@ contains !-------------------------------------------------------------------------------------------------- module subroutine plastic_none_init - integer :: & - Ninstance, & - p, & - NipcMyPhase - - write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_NONE_label//' init -+>>>' - - Ninstance = count(phase_plasticity == PLASTICITY_NONE_ID) - if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & - write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - - do p = 1, size(phase_plasticity) - if (phase_plasticity(p) /= PLASTICITY_NONE_ID) cycle - - NipcMyPhase = count(material_phaseAt == p) * discretization_nIP - call material_allocatePlasticState(p,NipcMyPhase,0,0,0) - - enddo + integer :: & + Ninstance, & + p, & + NipcMyPhase + + write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_NONE_label//' init -+>>>' + + Ninstance = count(phase_plasticity == PLASTICITY_NONE_ID) + if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & + write(6,'(a16,1x,i5,/)') '# instances:',Ninstance + + do p = 1, size(phase_plasticity) + if (phase_plasticity(p) /= PLASTICITY_NONE_ID) cycle + + NipcMyPhase = count(material_phaseAt == p) * discretization_nIP + call material_allocatePlasticState(p,NipcMyPhase,0,0,0) + + enddo end subroutine plastic_none_init diff --git a/src/constitutive_plastic_phenopowerlaw.f90 b/src/constitutive_plastic_phenopowerlaw.f90 index 3a330ef2a..63ff01c02 100644 --- a/src/constitutive_plastic_phenopowerlaw.f90 +++ b/src/constitutive_plastic_phenopowerlaw.f90 @@ -6,78 +6,78 @@ !-------------------------------------------------------------------------------------------------- submodule(constitutive) plastic_phenopowerlaw - enum, bind(c) - enumerator :: & - undefined_ID, & - resistance_slip_ID, & - accumulatedshear_slip_ID, & - shearrate_slip_ID, & - resolvedstress_slip_ID, & - resistance_twin_ID, & - accumulatedshear_twin_ID, & - shearrate_twin_ID, & - resolvedstress_twin_ID - end enum - - type :: tParameters - real(pReal) :: & - gdot0_slip, & !< reference shear strain rate for slip - gdot0_twin, & !< reference shear strain rate for twin - n_slip, & !< stress exponent for slip - n_twin, & !< stress exponent for twin - spr, & !< push-up factor for slip saturation due to twinning - twinB, & - twinC, & - twinD, & - twinE, & - h0_SlipSlip, & !< reference hardening slip - slip - h0_TwinSlip, & !< reference hardening twin - slip - h0_TwinTwin, & !< reference hardening twin - twin - a_slip, & - aTolResistance, & !< absolute tolerance for integration of xi - aTolShear, & !< absolute tolerance for integration of gamma - aTolTwinfrac !< absolute tolerance for integration of f - real(pReal), allocatable, dimension(:) :: & - xi_slip_0, & !< initial critical shear stress for slip - xi_twin_0, & !< initial critical shear stress for twin - xi_slip_sat, & !< maximum critical shear stress for slip - nonSchmidCoeff, & - H_int, & !< per family hardening activity (optional) - gamma_twin_char !< characteristic shear for twins - real(pReal), allocatable, dimension(:,:) :: & - interaction_SlipSlip, & !< slip resistance from slip activity - interaction_SlipTwin, & !< slip resistance from twin activity - interaction_TwinSlip, & !< twin resistance from slip activity - interaction_TwinTwin !< twin resistance from twin activity - real(pReal), allocatable, dimension(:,:,:) :: & - Schmid_slip, & - Schmid_twin, & - nonSchmid_pos, & - nonSchmid_neg - integer :: & - totalNslip, & !< total number of active slip system - totalNtwin !< total number of active twin systems - integer, allocatable, dimension(:) :: & - Nslip, & !< number of active slip systems for each family - Ntwin !< number of active twin systems for each family - integer(kind(undefined_ID)), allocatable, dimension(:) :: & - outputID !< ID of each post result output - end type tParameters - - type :: tPhenopowerlawState - real(pReal), pointer, dimension(:,:) :: & - xi_slip, & - xi_twin, & - gamma_slip, & - gamma_twin - end type tPhenopowerlawState - + enum, bind(c) + enumerator :: & + undefined_ID, & + resistance_slip_ID, & + accumulatedshear_slip_ID, & + shearrate_slip_ID, & + resolvedstress_slip_ID, & + resistance_twin_ID, & + accumulatedshear_twin_ID, & + shearrate_twin_ID, & + resolvedstress_twin_ID + end enum + + type :: tParameters + real(pReal) :: & + gdot0_slip, & !< reference shear strain rate for slip + gdot0_twin, & !< reference shear strain rate for twin + n_slip, & !< stress exponent for slip + n_twin, & !< stress exponent for twin + spr, & !< push-up factor for slip saturation due to twinning + twinB, & + twinC, & + twinD, & + twinE, & + h0_SlipSlip, & !< reference hardening slip - slip + h0_TwinSlip, & !< reference hardening twin - slip + h0_TwinTwin, & !< reference hardening twin - twin + a_slip, & + aTolResistance, & !< absolute tolerance for integration of xi + aTolShear, & !< absolute tolerance for integration of gamma + aTolTwinfrac !< absolute tolerance for integration of f + real(pReal), allocatable, dimension(:) :: & + xi_slip_0, & !< initial critical shear stress for slip + xi_twin_0, & !< initial critical shear stress for twin + xi_slip_sat, & !< maximum critical shear stress for slip + nonSchmidCoeff, & + H_int, & !< per family hardening activity (optional) + gamma_twin_char !< characteristic shear for twins + real(pReal), allocatable, dimension(:,:) :: & + interaction_SlipSlip, & !< slip resistance from slip activity + interaction_SlipTwin, & !< slip resistance from twin activity + interaction_TwinSlip, & !< twin resistance from slip activity + interaction_TwinTwin !< twin resistance from twin activity + real(pReal), allocatable, dimension(:,:,:) :: & + Schmid_slip, & + Schmid_twin, & + nonSchmid_pos, & + nonSchmid_neg + integer :: & + totalNslip, & !< total number of active slip system + totalNtwin !< total number of active twin systems + integer, allocatable, dimension(:) :: & + Nslip, & !< number of active slip systems for each family + Ntwin !< number of active twin systems for each family + integer(kind(undefined_ID)), allocatable, dimension(:) :: & + outputID !< ID of each post result output + end type tParameters + + type :: tPhenopowerlawState + real(pReal), pointer, dimension(:,:) :: & + xi_slip, & + xi_twin, & + gamma_slip, & + gamma_twin + end type tPhenopowerlawState + !-------------------------------------------------------------------------------------------------- ! containers for parameters and state - type(tParameters), allocatable, dimension(:) :: param - type(tPhenopowerlawState), allocatable, dimension(:) :: & - dotState, & - state + type(tParameters), allocatable, dimension(:) :: param + type(tPhenopowerlawState), allocatable, dimension(:) :: & + dotState, & + state contains @@ -88,242 +88,242 @@ contains !-------------------------------------------------------------------------------------------------- module subroutine plastic_phenopowerlaw_init - integer :: & - Ninstance, & - p, i, & - NipcMyPhase, outputSize, & - sizeState, sizeDotState, & - startIndex, endIndex - - integer(kind(undefined_ID)) :: & - outputID - - character(len=pStringLen) :: & - extmsg = '' - character(len=pStringLen), dimension(:), allocatable :: & - outputs - - write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_PHENOPOWERLAW_label//' init -+>>>' - - Ninstance = count(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID) - if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & - write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - - allocate(param(Ninstance)) - allocate(state(Ninstance)) - allocate(dotState(Ninstance)) - - do p = 1, size(phase_plasticity) - if (phase_plasticity(p) /= PLASTICITY_PHENOPOWERLAW_ID) cycle - associate(prm => param(phase_plasticityInstance(p)), & - dot => dotState(phase_plasticityInstance(p)), & - stt => state(phase_plasticityInstance(p)), & - config => config_phase(p)) - + integer :: & + Ninstance, & + p, i, & + NipcMyPhase, outputSize, & + sizeState, sizeDotState, & + startIndex, endIndex + + integer(kind(undefined_ID)) :: & + outputID + + character(len=pStringLen) :: & + extmsg = '' + character(len=pStringLen), dimension(:), allocatable :: & + outputs + + write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_PHENOPOWERLAW_label//' init -+>>>' + + Ninstance = count(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID) + if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & + write(6,'(a16,1x,i5,/)') '# instances:',Ninstance + + allocate(param(Ninstance)) + allocate(state(Ninstance)) + allocate(dotState(Ninstance)) + + do p = 1, size(phase_plasticity) + if (phase_plasticity(p) /= PLASTICITY_PHENOPOWERLAW_ID) cycle + associate(prm => param(phase_plasticityInstance(p)), & + dot => dotState(phase_plasticityInstance(p)), & + stt => state(phase_plasticityInstance(p)), & + config => config_phase(p)) + !-------------------------------------------------------------------------------------------------- ! optional parameters that need to be defined - prm%twinB = config%getFloat('twin_b',defaultVal=1.0_pReal) - prm%twinC = config%getFloat('twin_c',defaultVal=0.0_pReal) - prm%twinD = config%getFloat('twin_d',defaultVal=0.0_pReal) - prm%twinE = config%getFloat('twin_e',defaultVal=0.0_pReal) - - prm%aTolResistance = config%getFloat('atol_resistance',defaultVal=1.0_pReal) - prm%aTolShear = config%getFloat('atol_shear', defaultVal=1.0e-6_pReal) - prm%aTolTwinfrac = config%getFloat('atol_twinfrac', defaultVal=1.0e-6_pReal) - - ! sanity checks - if (prm%aTolResistance <= 0.0_pReal) extmsg = trim(extmsg)//' aTolresistance' - if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//' aTolShear' - if (prm%aTolTwinfrac <= 0.0_pReal) extmsg = trim(extmsg)//' atoltwinfrac' - + prm%twinB = config%getFloat('twin_b',defaultVal=1.0_pReal) + prm%twinC = config%getFloat('twin_c',defaultVal=0.0_pReal) + prm%twinD = config%getFloat('twin_d',defaultVal=0.0_pReal) + prm%twinE = config%getFloat('twin_e',defaultVal=0.0_pReal) + + prm%aTolResistance = config%getFloat('atol_resistance',defaultVal=1.0_pReal) + prm%aTolShear = config%getFloat('atol_shear', defaultVal=1.0e-6_pReal) + prm%aTolTwinfrac = config%getFloat('atol_twinfrac', defaultVal=1.0e-6_pReal) + + ! sanity checks + if (prm%aTolResistance <= 0.0_pReal) extmsg = trim(extmsg)//' aTolresistance' + if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//' aTolShear' + if (prm%aTolTwinfrac <= 0.0_pReal) extmsg = trim(extmsg)//' atoltwinfrac' + !-------------------------------------------------------------------------------------------------- ! slip related parameters - prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) - prm%totalNslip = sum(prm%Nslip) - slipActive: if (prm%totalNslip > 0) then - prm%Schmid_slip = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),& - config%getFloat('c/a',defaultVal=0.0_pReal)) - - if(trim(config%getString('lattice_structure')) == 'bcc') then - prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',& - defaultVal = emptyRealArray) - prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1) - prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1) - else - allocate(prm%nonSchmidCoeff(0)) - prm%nonSchmid_pos = prm%Schmid_slip - prm%nonSchmid_neg = prm%Schmid_slip - endif - prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(prm%Nslip, & - config%getFloats('interaction_slipslip'), & - config%getString('lattice_structure')) - - prm%xi_slip_0 = config%getFloats('tau0_slip', requiredSize=size(prm%Nslip)) - prm%xi_slip_sat = config%getFloats('tausat_slip', requiredSize=size(prm%Nslip)) - prm%H_int = config%getFloats('h_int', requiredSize=size(prm%Nslip), & - defaultVal=[(0.0_pReal,i=1,size(prm%Nslip))]) - - prm%gdot0_slip = config%getFloat('gdot0_slip') - prm%n_slip = config%getFloat('n_slip') - prm%a_slip = config%getFloat('a_slip') - prm%h0_SlipSlip = config%getFloat('h0_slipslip') - - ! expand: family => system - prm%xi_slip_0 = math_expand(prm%xi_slip_0, prm%Nslip) - prm%xi_slip_sat = math_expand(prm%xi_slip_sat,prm%Nslip) - prm%H_int = math_expand(prm%H_int, prm%Nslip) - - ! sanity checks - if ( prm%gdot0_slip <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0_slip' - if ( prm%a_slip <= 0.0_pReal) extmsg = trim(extmsg)//' a_slip' - if ( prm%n_slip <= 0.0_pReal) extmsg = trim(extmsg)//' n_slip' - if (any(prm%xi_slip_0 <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_slip_0' - if (any(prm%xi_slip_sat <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_slip_sat' - else slipActive - allocate(prm%interaction_SlipSlip(0,0)) - allocate(prm%xi_slip_0(0)) - endif slipActive - + prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) + prm%totalNslip = sum(prm%Nslip) + slipActive: if (prm%totalNslip > 0) then + prm%Schmid_slip = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),& + config%getFloat('c/a',defaultVal=0.0_pReal)) + + if(trim(config%getString('lattice_structure')) == 'bcc') then + prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',& + defaultVal = emptyRealArray) + prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1) + prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1) + else + allocate(prm%nonSchmidCoeff(0)) + prm%nonSchmid_pos = prm%Schmid_slip + prm%nonSchmid_neg = prm%Schmid_slip + endif + prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(prm%Nslip, & + config%getFloats('interaction_slipslip'), & + config%getString('lattice_structure')) + + prm%xi_slip_0 = config%getFloats('tau0_slip', requiredSize=size(prm%Nslip)) + prm%xi_slip_sat = config%getFloats('tausat_slip', requiredSize=size(prm%Nslip)) + prm%H_int = config%getFloats('h_int', requiredSize=size(prm%Nslip), & + defaultVal=[(0.0_pReal,i=1,size(prm%Nslip))]) + + prm%gdot0_slip = config%getFloat('gdot0_slip') + prm%n_slip = config%getFloat('n_slip') + prm%a_slip = config%getFloat('a_slip') + prm%h0_SlipSlip = config%getFloat('h0_slipslip') + + ! expand: family => system + prm%xi_slip_0 = math_expand(prm%xi_slip_0, prm%Nslip) + prm%xi_slip_sat = math_expand(prm%xi_slip_sat,prm%Nslip) + prm%H_int = math_expand(prm%H_int, prm%Nslip) + + ! sanity checks + if ( prm%gdot0_slip <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0_slip' + if ( prm%a_slip <= 0.0_pReal) extmsg = trim(extmsg)//' a_slip' + if ( prm%n_slip <= 0.0_pReal) extmsg = trim(extmsg)//' n_slip' + if (any(prm%xi_slip_0 <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_slip_0' + if (any(prm%xi_slip_sat <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_slip_sat' + else slipActive + allocate(prm%interaction_SlipSlip(0,0)) + allocate(prm%xi_slip_0(0)) + endif slipActive + !-------------------------------------------------------------------------------------------------- ! twin related parameters - prm%Ntwin = config%getInts('ntwin', defaultVal=emptyIntArray) - prm%totalNtwin = sum(prm%Ntwin) - twinActive: if (prm%totalNtwin > 0) then - prm%Schmid_twin = lattice_SchmidMatrix_twin(prm%Ntwin,config%getString('lattice_structure'),& - config%getFloat('c/a',defaultVal=0.0_pReal)) - prm%interaction_TwinTwin = lattice_interaction_TwinByTwin(prm%Ntwin,& - config%getFloats('interaction_twintwin'), & - config%getString('lattice_structure')) - prm%gamma_twin_char = lattice_characteristicShear_twin(prm%Ntwin,config%getString('lattice_structure'),& - config%getFloat('c/a')) - - prm%xi_twin_0 = config%getFloats('tau0_twin',requiredSize=size(prm%Ntwin)) - - prm%gdot0_twin = config%getFloat('gdot0_twin') - prm%n_twin = config%getFloat('n_twin') - prm%spr = config%getFloat('s_pr') - prm%h0_TwinTwin = config%getFloat('h0_twintwin') - - ! expand: family => system - prm%xi_twin_0 = math_expand(prm%xi_twin_0, prm%Ntwin) - - ! sanity checks - if (prm%gdot0_twin <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0_twin' - if (prm%n_twin <= 0.0_pReal) extmsg = trim(extmsg)//' n_twin' - else twinActive - allocate(prm%interaction_TwinTwin(0,0)) - allocate(prm%xi_twin_0(0)) - allocate(prm%gamma_twin_char(0)) - endif twinActive - + prm%Ntwin = config%getInts('ntwin', defaultVal=emptyIntArray) + prm%totalNtwin = sum(prm%Ntwin) + twinActive: if (prm%totalNtwin > 0) then + prm%Schmid_twin = lattice_SchmidMatrix_twin(prm%Ntwin,config%getString('lattice_structure'),& + config%getFloat('c/a',defaultVal=0.0_pReal)) + prm%interaction_TwinTwin = lattice_interaction_TwinByTwin(prm%Ntwin,& + config%getFloats('interaction_twintwin'), & + config%getString('lattice_structure')) + prm%gamma_twin_char = lattice_characteristicShear_twin(prm%Ntwin,config%getString('lattice_structure'),& + config%getFloat('c/a')) + + prm%xi_twin_0 = config%getFloats('tau0_twin',requiredSize=size(prm%Ntwin)) + + prm%gdot0_twin = config%getFloat('gdot0_twin') + prm%n_twin = config%getFloat('n_twin') + prm%spr = config%getFloat('s_pr') + prm%h0_TwinTwin = config%getFloat('h0_twintwin') + + ! expand: family => system + prm%xi_twin_0 = math_expand(prm%xi_twin_0, prm%Ntwin) + + ! sanity checks + if (prm%gdot0_twin <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0_twin' + if (prm%n_twin <= 0.0_pReal) extmsg = trim(extmsg)//' n_twin' + else twinActive + allocate(prm%interaction_TwinTwin(0,0)) + allocate(prm%xi_twin_0(0)) + allocate(prm%gamma_twin_char(0)) + endif twinActive + !-------------------------------------------------------------------------------------------------- ! slip-twin related parameters - slipAndTwinActive: if (prm%totalNslip > 0 .and. prm%totalNtwin > 0) then - prm%h0_TwinSlip = config%getFloat('h0_twinslip') - prm%interaction_SlipTwin = lattice_interaction_SlipByTwin(prm%Nslip,prm%Ntwin,& - config%getFloats('interaction_sliptwin'), & - config%getString('lattice_structure')) - prm%interaction_TwinSlip = lattice_interaction_TwinBySlip(prm%Ntwin,prm%Nslip,& - config%getFloats('interaction_twinslip'), & - config%getString('lattice_structure')) - else slipAndTwinActive - allocate(prm%interaction_SlipTwin(prm%TotalNslip,prm%TotalNtwin)) ! at least one dimension is 0 - allocate(prm%interaction_TwinSlip(prm%TotalNtwin,prm%TotalNslip)) ! at least one dimension is 0 - prm%h0_TwinSlip = 0.0_pReal - endif slipAndTwinActive - + slipAndTwinActive: if (prm%totalNslip > 0 .and. prm%totalNtwin > 0) then + prm%h0_TwinSlip = config%getFloat('h0_twinslip') + prm%interaction_SlipTwin = lattice_interaction_SlipByTwin(prm%Nslip,prm%Ntwin,& + config%getFloats('interaction_sliptwin'), & + config%getString('lattice_structure')) + prm%interaction_TwinSlip = lattice_interaction_TwinBySlip(prm%Ntwin,prm%Nslip,& + config%getFloats('interaction_twinslip'), & + config%getString('lattice_structure')) + else slipAndTwinActive + allocate(prm%interaction_SlipTwin(prm%TotalNslip,prm%TotalNtwin)) ! at least one dimension is 0 + allocate(prm%interaction_TwinSlip(prm%TotalNtwin,prm%TotalNslip)) ! at least one dimension is 0 + prm%h0_TwinSlip = 0.0_pReal + endif slipAndTwinActive + !-------------------------------------------------------------------------------------------------- ! exit if any parameter is out of range - if (extmsg /= '') & - call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_PHENOPOWERLAW_label//')') - + if (extmsg /= '') & + call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_PHENOPOWERLAW_label//')') + !-------------------------------------------------------------------------------------------------- ! output pararameters - outputs = config%getStrings('(output)',defaultVal=emptyStringArray) - allocate(prm%outputID(0)) - do i=1, size(outputs) - outputID = undefined_ID - select case(outputs(i)) - - case ('resistance_slip') - outputID = merge(resistance_slip_ID,undefined_ID,prm%totalNslip>0) - outputSize = prm%totalNslip - case ('accumulatedshear_slip') - outputID = merge(accumulatedshear_slip_ID,undefined_ID,prm%totalNslip>0) - outputSize = prm%totalNslip - case ('shearrate_slip') - outputID = merge(shearrate_slip_ID,undefined_ID,prm%totalNslip>0) - outputSize = prm%totalNslip - case ('resolvedstress_slip') - outputID = merge(resolvedstress_slip_ID,undefined_ID,prm%totalNslip>0) - outputSize = prm%totalNslip - - case ('resistance_twin') - outputID = merge(resistance_twin_ID,undefined_ID,prm%totalNtwin>0) - outputSize = prm%totalNtwin - case ('accumulatedshear_twin') - outputID = merge(accumulatedshear_twin_ID,undefined_ID,prm%totalNtwin>0) - outputSize = prm%totalNtwin - case ('shearrate_twin') - outputID = merge(shearrate_twin_ID,undefined_ID,prm%totalNtwin>0) - outputSize = prm%totalNtwin - case ('resolvedstress_twin') - outputID = merge(resolvedstress_twin_ID,undefined_ID,prm%totalNtwin>0) - outputSize = prm%totalNtwin - - end select - - if (outputID /= undefined_ID) then - prm%outputID = [prm%outputID, outputID] - endif - - enddo - + outputs = config%getStrings('(output)',defaultVal=emptyStringArray) + allocate(prm%outputID(0)) + do i=1, size(outputs) + outputID = undefined_ID + select case(outputs(i)) + + case ('resistance_slip') + outputID = merge(resistance_slip_ID,undefined_ID,prm%totalNslip>0) + outputSize = prm%totalNslip + case ('accumulatedshear_slip') + outputID = merge(accumulatedshear_slip_ID,undefined_ID,prm%totalNslip>0) + outputSize = prm%totalNslip + case ('shearrate_slip') + outputID = merge(shearrate_slip_ID,undefined_ID,prm%totalNslip>0) + outputSize = prm%totalNslip + case ('resolvedstress_slip') + outputID = merge(resolvedstress_slip_ID,undefined_ID,prm%totalNslip>0) + outputSize = prm%totalNslip + + case ('resistance_twin') + outputID = merge(resistance_twin_ID,undefined_ID,prm%totalNtwin>0) + outputSize = prm%totalNtwin + case ('accumulatedshear_twin') + outputID = merge(accumulatedshear_twin_ID,undefined_ID,prm%totalNtwin>0) + outputSize = prm%totalNtwin + case ('shearrate_twin') + outputID = merge(shearrate_twin_ID,undefined_ID,prm%totalNtwin>0) + outputSize = prm%totalNtwin + case ('resolvedstress_twin') + outputID = merge(resolvedstress_twin_ID,undefined_ID,prm%totalNtwin>0) + outputSize = prm%totalNtwin + + end select + + if (outputID /= undefined_ID) then + prm%outputID = [prm%outputID, outputID] + endif + + enddo + !-------------------------------------------------------------------------------------------------- ! allocate state arrays - NipcMyPhase = count(material_phaseAt == p) * discretization_nIP - sizeDotState = size(['tau_slip ','gamma_slip']) * prm%totalNslip & - + size(['tau_twin ','gamma_twin']) * prm%totalNtwin - sizeState = sizeDotState - - call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0) - + NipcMyPhase = count(material_phaseAt == p) * discretization_nIP + sizeDotState = size(['tau_slip ','gamma_slip']) * prm%totalNslip & + + size(['tau_twin ','gamma_twin']) * prm%totalNtwin + sizeState = sizeDotState + + call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0) + !-------------------------------------------------------------------------------------------------- ! locally defined state aliases and initialization of state0 and aTolState - startIndex = 1 - endIndex = prm%totalNslip - stt%xi_slip => plasticState(p)%state (startIndex:endIndex,:) - stt%xi_slip = spread(prm%xi_slip_0, 2, NipcMyPhase) - dot%xi_slip => plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance - - startIndex = endIndex + 1 - endIndex = endIndex + prm%totalNtwin - stt%xi_twin => plasticState(p)%state (startIndex:endIndex,:) - stt%xi_twin = spread(prm%xi_twin_0, 2, NipcMyPhase) - dot%xi_twin => plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance - - startIndex = endIndex + 1 - endIndex = endIndex + prm%totalNslip - stt%gamma_slip => plasticState(p)%state (startIndex:endIndex,:) - dot%gamma_slip => plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear - ! global alias - plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%accumulatedSlip => plasticState(p)%state(startIndex:endIndex,:) - - startIndex = endIndex + 1 - endIndex = endIndex + prm%totalNtwin - stt%gamma_twin => plasticState(p)%state (startIndex:endIndex,:) - dot%gamma_twin => plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear - - plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally - - end associate - - enddo + startIndex = 1 + endIndex = prm%totalNslip + stt%xi_slip => plasticState(p)%state (startIndex:endIndex,:) + stt%xi_slip = spread(prm%xi_slip_0, 2, NipcMyPhase) + dot%xi_slip => plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance + + startIndex = endIndex + 1 + endIndex = endIndex + prm%totalNtwin + stt%xi_twin => plasticState(p)%state (startIndex:endIndex,:) + stt%xi_twin = spread(prm%xi_twin_0, 2, NipcMyPhase) + dot%xi_twin => plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance + + startIndex = endIndex + 1 + endIndex = endIndex + prm%totalNslip + stt%gamma_slip => plasticState(p)%state (startIndex:endIndex,:) + dot%gamma_slip => plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear + ! global alias + plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%accumulatedSlip => plasticState(p)%state(startIndex:endIndex,:) + + startIndex = endIndex + 1 + endIndex = endIndex + prm%totalNtwin + stt%gamma_twin => plasticState(p)%state (startIndex:endIndex,:) + dot%gamma_twin => plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear + + plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally + + end associate + + enddo end subroutine plastic_phenopowerlaw_init @@ -335,48 +335,48 @@ end subroutine plastic_phenopowerlaw_init !-------------------------------------------------------------------------------------------------- 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 - - integer :: & - i,k,l,m,n - real(pReal), dimension(param(instance)%totalNslip) :: & - gdot_slip_pos,gdot_slip_neg, & - dgdot_dtauslip_pos,dgdot_dtauslip_neg - real(pReal), dimension(param(instance)%totalNtwin) :: & - gdot_twin,dgdot_dtautwin - - Lp = 0.0_pReal - dLp_dMp = 0.0_pReal - - associate(prm => param(instance)) - - call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg) - slipSystems: do i = 1, prm%totalNslip - Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%Schmid_slip(1:3,1:3,i) - forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & - + dgdot_dtauslip_pos(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_pos(m,n,i) & - + dgdot_dtauslip_neg(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_neg(m,n,i) - enddo slipSystems - - call kinetics_twin(Mp,instance,of,gdot_twin,dgdot_dtautwin) - twinSystems: do i = 1, prm%totalNtwin - Lp = Lp + gdot_twin(i)*prm%Schmid_twin(1:3,1:3,i) - forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & - + dgdot_dtautwin(i)*prm%Schmid_twin(k,l,i)*prm%Schmid_twin(m,n,i) - enddo twinSystems - - end associate + 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 + + integer :: & + i,k,l,m,n + real(pReal), dimension(param(instance)%totalNslip) :: & + gdot_slip_pos,gdot_slip_neg, & + dgdot_dtauslip_pos,dgdot_dtauslip_neg + real(pReal), dimension(param(instance)%totalNtwin) :: & + gdot_twin,dgdot_dtautwin + + Lp = 0.0_pReal + dLp_dMp = 0.0_pReal + + associate(prm => param(instance)) + + call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg) + slipSystems: do i = 1, prm%totalNslip + Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%Schmid_slip(1:3,1:3,i) + forall (k=1:3,l=1:3,m=1:3,n=1:3) & + dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + + dgdot_dtauslip_pos(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_pos(m,n,i) & + + dgdot_dtauslip_neg(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_neg(m,n,i) + enddo slipSystems + + call kinetics_twin(Mp,instance,of,gdot_twin,dgdot_dtautwin) + twinSystems: do i = 1, prm%totalNtwin + Lp = Lp + gdot_twin(i)*prm%Schmid_twin(1:3,1:3,i) + forall (k=1:3,l=1:3,m=1:3,n=1:3) & + dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + + dgdot_dtautwin(i)*prm%Schmid_twin(k,l,i)*prm%Schmid_twin(m,n,i) + enddo twinSystems + + end associate end subroutine plastic_phenopowerlaw_LpAndItsTangent @@ -386,53 +386,53 @@ end subroutine plastic_phenopowerlaw_LpAndItsTangent !-------------------------------------------------------------------------------------------------- module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - - real(pReal) :: & - c_SlipSlip,c_TwinSlip,c_TwinTwin, & - xi_slip_sat_offset,& - sumGamma,sumF - real(pReal), dimension(param(instance)%totalNslip) :: & - left_SlipSlip,right_SlipSlip, & - gdot_slip_pos,gdot_slip_neg - - associate(prm => param(instance), stt => state(instance), dot => dotState(instance)) - - sumGamma = sum(stt%gamma_slip(:,of)) - sumF = sum(stt%gamma_twin(:,of)/prm%gamma_twin_char) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + + real(pReal) :: & + c_SlipSlip,c_TwinSlip,c_TwinTwin, & + xi_slip_sat_offset,& + sumGamma,sumF + real(pReal), dimension(param(instance)%totalNslip) :: & + left_SlipSlip,right_SlipSlip, & + gdot_slip_pos,gdot_slip_neg + + associate(prm => param(instance), stt => state(instance), dot => dotState(instance)) + + sumGamma = sum(stt%gamma_slip(:,of)) + sumF = sum(stt%gamma_twin(:,of)/prm%gamma_twin_char) !-------------------------------------------------------------------------------------------------- ! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices - c_SlipSlip = prm%h0_slipslip * (1.0_pReal + prm%twinC*sumF** prm%twinB) - c_TwinSlip = prm%h0_TwinSlip * sumGamma**prm%twinE - c_TwinTwin = prm%h0_TwinTwin * sumF**prm%twinD - + c_SlipSlip = prm%h0_slipslip * (1.0_pReal + prm%twinC*sumF** prm%twinB) + c_TwinSlip = prm%h0_TwinSlip * sumGamma**prm%twinE + c_TwinTwin = prm%h0_TwinTwin * sumF**prm%twinD + !-------------------------------------------------------------------------------------------------- ! calculate left and right vectors - left_SlipSlip = 1.0_pReal + prm%H_int - xi_slip_sat_offset = prm%spr*sqrt(sumF) - right_SlipSlip = abs(1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset)) **prm%a_slip & - * sign(1.0_pReal,1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset)) - + left_SlipSlip = 1.0_pReal + prm%H_int + xi_slip_sat_offset = prm%spr*sqrt(sumF) + right_SlipSlip = abs(1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset)) **prm%a_slip & + * sign(1.0_pReal,1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+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,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)) + !-------------------------------------------------------------------------------------------------- ! hardening - dot%xi_slip(:,of) = c_SlipSlip * left_SlipSlip * & - matmul(prm%interaction_SlipSlip,dot%gamma_slip(:,of)*right_SlipSlip) & - + matmul(prm%interaction_SlipTwin,dot%gamma_twin(:,of)) - - dot%xi_twin(:,of) = c_TwinSlip * matmul(prm%interaction_TwinSlip,dot%gamma_slip(:,of)) & - + c_TwinTwin * matmul(prm%interaction_TwinTwin,dot%gamma_twin(:,of)) - end associate + dot%xi_slip(:,of) = c_SlipSlip * left_SlipSlip * & + matmul(prm%interaction_SlipSlip,dot%gamma_slip(:,of)*right_SlipSlip) & + + matmul(prm%interaction_SlipTwin,dot%gamma_twin(:,of)) + + dot%xi_twin(:,of) = c_TwinSlip * matmul(prm%interaction_TwinSlip,dot%gamma_slip(:,of)) & + + c_TwinTwin * matmul(prm%interaction_TwinTwin,dot%gamma_twin(:,of)) + end associate end subroutine plastic_phenopowerlaw_dotState @@ -481,64 +481,64 @@ end subroutine plastic_phenopowerlaw_results pure subroutine kinetics_slip(Mp,instance,of, & 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 - - real(pReal), intent(out), dimension(param(instance)%totalNslip) :: & - gdot_slip_pos, & - gdot_slip_neg - real(pReal), intent(out), optional, dimension(param(instance)%totalNslip) :: & - dgdot_dtau_slip_pos, & - dgdot_dtau_slip_neg - - real(pReal), dimension(param(instance)%totalNslip) :: & - tau_slip_pos, & - tau_slip_neg - integer :: i - logical :: nonSchmidActive - - associate(prm => param(instance), stt => state(instance)) - - nonSchmidActive = size(prm%nonSchmidCoeff) > 0 - - do i = 1, prm%totalNslip - tau_slip_pos(i) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i)) - tau_slip_neg(i) = merge(math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i)), & - 0.0_pReal, nonSchmidActive) - enddo - - where(dNeq0(tau_slip_pos)) - gdot_slip_pos = prm%gdot0_slip * merge(0.5_pReal,1.0_pReal, nonSchmidActive) & ! 1/2 if non-Schmid active - * sign(abs(tau_slip_pos/stt%xi_slip(:,of))**prm%n_slip, tau_slip_pos) - else where - gdot_slip_pos = 0.0_pReal - end where - - where(dNeq0(tau_slip_neg)) - gdot_slip_neg = prm%gdot0_slip * 0.5_pReal & ! only used if non-Schmid active, always 1/2 - * sign(abs(tau_slip_neg/stt%xi_slip(:,of))**prm%n_slip, tau_slip_neg) - else where - gdot_slip_neg = 0.0_pReal - end where - - if (present(dgdot_dtau_slip_pos)) then - where(dNeq0(gdot_slip_pos)) - dgdot_dtau_slip_pos = gdot_slip_pos*prm%n_slip/tau_slip_pos - else where - dgdot_dtau_slip_pos = 0.0_pReal - end where - endif - if (present(dgdot_dtau_slip_neg)) then - where(dNeq0(gdot_slip_neg)) - dgdot_dtau_slip_neg = gdot_slip_neg*prm%n_slip/tau_slip_neg - else where - dgdot_dtau_slip_neg = 0.0_pReal - end where - endif - end associate + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + + real(pReal), intent(out), dimension(param(instance)%totalNslip) :: & + gdot_slip_pos, & + gdot_slip_neg + real(pReal), intent(out), optional, dimension(param(instance)%totalNslip) :: & + dgdot_dtau_slip_pos, & + dgdot_dtau_slip_neg + + real(pReal), dimension(param(instance)%totalNslip) :: & + tau_slip_pos, & + tau_slip_neg + integer :: i + logical :: nonSchmidActive + + associate(prm => param(instance), stt => state(instance)) + + nonSchmidActive = size(prm%nonSchmidCoeff) > 0 + + do i = 1, prm%totalNslip + tau_slip_pos(i) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i)) + tau_slip_neg(i) = merge(math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i)), & + 0.0_pReal, nonSchmidActive) + enddo + + where(dNeq0(tau_slip_pos)) + gdot_slip_pos = prm%gdot0_slip * merge(0.5_pReal,1.0_pReal, nonSchmidActive) & ! 1/2 if non-Schmid active + * sign(abs(tau_slip_pos/stt%xi_slip(:,of))**prm%n_slip, tau_slip_pos) + else where + gdot_slip_pos = 0.0_pReal + end where + + where(dNeq0(tau_slip_neg)) + gdot_slip_neg = prm%gdot0_slip * 0.5_pReal & ! only used if non-Schmid active, always 1/2 + * sign(abs(tau_slip_neg/stt%xi_slip(:,of))**prm%n_slip, tau_slip_neg) + else where + gdot_slip_neg = 0.0_pReal + end where + + if (present(dgdot_dtau_slip_pos)) then + where(dNeq0(gdot_slip_pos)) + dgdot_dtau_slip_pos = gdot_slip_pos*prm%n_slip/tau_slip_pos + else where + dgdot_dtau_slip_pos = 0.0_pReal + end where + endif + if (present(dgdot_dtau_slip_neg)) then + where(dNeq0(gdot_slip_neg)) + dgdot_dtau_slip_neg = gdot_slip_neg*prm%n_slip/tau_slip_neg + else where + dgdot_dtau_slip_neg = 0.0_pReal + end where + endif + end associate end subroutine kinetics_slip @@ -553,43 +553,43 @@ end subroutine kinetics_slip pure subroutine kinetics_twin(Mp,instance,of,& gdot_twin,dgdot_dtau_twin) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - - real(pReal), dimension(param(instance)%totalNtwin), intent(out) :: & - gdot_twin - real(pReal), dimension(param(instance)%totalNtwin), intent(out), optional :: & - dgdot_dtau_twin - - real(pReal), dimension(param(instance)%totalNtwin) :: & - tau_twin - integer :: i - - associate(prm => param(instance), stt => state(instance)) - - do i = 1, prm%totalNtwin - tau_twin(i) = math_mul33xx33(Mp,prm%Schmid_twin(1:3,1:3,i)) - enddo - - where(tau_twin > 0.0_pReal) - gdot_twin = (1.0_pReal-sum(stt%gamma_twin(:,of)/prm%gamma_twin_char)) & ! only twin in untwinned volume fraction - * prm%gdot0_twin*(abs(tau_twin)/stt%xi_twin(:,of))**prm%n_twin - else where - gdot_twin = 0.0_pReal - end where - - if (present(dgdot_dtau_twin)) then - where(dNeq0(gdot_twin)) - dgdot_dtau_twin = gdot_twin*prm%n_twin/tau_twin - else where - dgdot_dtau_twin = 0.0_pReal - end where - endif - - end associate + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + + real(pReal), dimension(param(instance)%totalNtwin), intent(out) :: & + gdot_twin + real(pReal), dimension(param(instance)%totalNtwin), intent(out), optional :: & + dgdot_dtau_twin + + real(pReal), dimension(param(instance)%totalNtwin) :: & + tau_twin + integer :: i + + associate(prm => param(instance), stt => state(instance)) + + do i = 1, prm%totalNtwin + tau_twin(i) = math_mul33xx33(Mp,prm%Schmid_twin(1:3,1:3,i)) + enddo + + where(tau_twin > 0.0_pReal) + gdot_twin = (1.0_pReal-sum(stt%gamma_twin(:,of)/prm%gamma_twin_char)) & ! only twin in untwinned volume fraction + * prm%gdot0_twin*(abs(tau_twin)/stt%xi_twin(:,of))**prm%n_twin + else where + gdot_twin = 0.0_pReal + end where + + if (present(dgdot_dtau_twin)) then + where(dNeq0(gdot_twin)) + dgdot_dtau_twin = gdot_twin*prm%n_twin/tau_twin + else where + dgdot_dtau_twin = 0.0_pReal + end where + endif + + end associate end subroutine kinetics_twin From 042d09a73025a597d309d0457decfce47028b8d9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 31 Jan 2020 21:43:12 +0100 Subject: [PATCH 16/20] names from paper --- src/constitutive_plastic_phenopowerlaw.f90 | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/constitutive_plastic_phenopowerlaw.f90 b/src/constitutive_plastic_phenopowerlaw.f90 index 63ff01c02..e14cb830d 100644 --- a/src/constitutive_plastic_phenopowerlaw.f90 +++ b/src/constitutive_plastic_phenopowerlaw.f90 @@ -26,10 +26,10 @@ submodule(constitutive) plastic_phenopowerlaw n_slip, & !< stress exponent for slip n_twin, & !< stress exponent for twin spr, & !< push-up factor for slip saturation due to twinning - twinB, & - twinC, & - twinD, & - twinE, & + c_2, & + c_1, & + c_4, & + c_3, & h0_SlipSlip, & !< reference hardening slip - slip h0_TwinSlip, & !< reference hardening twin - slip h0_TwinTwin, & !< reference hardening twin - twin @@ -122,10 +122,10 @@ module subroutine plastic_phenopowerlaw_init !-------------------------------------------------------------------------------------------------- ! optional parameters that need to be defined - prm%twinB = config%getFloat('twin_b',defaultVal=1.0_pReal) - prm%twinC = config%getFloat('twin_c',defaultVal=0.0_pReal) - prm%twinD = config%getFloat('twin_d',defaultVal=0.0_pReal) - prm%twinE = config%getFloat('twin_e',defaultVal=0.0_pReal) + prm%c_2 = config%getFloat('twin_b',defaultVal=1.0_pReal) + prm%c_1 = config%getFloat('twin_c',defaultVal=0.0_pReal) + prm%c_4 = config%getFloat('twin_d',defaultVal=0.0_pReal) + prm%c_3 = config%getFloat('twin_e',defaultVal=0.0_pReal) prm%aTolResistance = config%getFloat('atol_resistance',defaultVal=1.0_pReal) prm%aTolShear = config%getFloat('atol_shear', defaultVal=1.0e-6_pReal) @@ -407,9 +407,9 @@ module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) !-------------------------------------------------------------------------------------------------- ! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices - c_SlipSlip = prm%h0_slipslip * (1.0_pReal + prm%twinC*sumF** prm%twinB) - c_TwinSlip = prm%h0_TwinSlip * sumGamma**prm%twinE - c_TwinTwin = prm%h0_TwinTwin * sumF**prm%twinD + c_SlipSlip = prm%h0_slipslip * (1.0_pReal + prm%c_1*sumF** prm%c_2) + c_TwinSlip = prm%h0_TwinSlip * sumGamma**prm%c_3 + c_TwinTwin = prm%h0_TwinTwin * sumF**prm%c_4 !-------------------------------------------------------------------------------------------------- ! calculate left and right vectors From 796fd9a7743baf5578cecdb0ba750c1c753c3e38 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 31 Jan 2020 21:43:45 +0100 Subject: [PATCH 17/20] natural order --- src/constitutive_plastic_phenopowerlaw.f90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/constitutive_plastic_phenopowerlaw.f90 b/src/constitutive_plastic_phenopowerlaw.f90 index e14cb830d..95fc2b46b 100644 --- a/src/constitutive_plastic_phenopowerlaw.f90 +++ b/src/constitutive_plastic_phenopowerlaw.f90 @@ -26,10 +26,10 @@ submodule(constitutive) plastic_phenopowerlaw n_slip, & !< stress exponent for slip n_twin, & !< stress exponent for twin spr, & !< push-up factor for slip saturation due to twinning - c_2, & c_1, & - c_4, & + c_2, & c_3, & + c_4, & h0_SlipSlip, & !< reference hardening slip - slip h0_TwinSlip, & !< reference hardening twin - slip h0_TwinTwin, & !< reference hardening twin - twin @@ -122,10 +122,10 @@ module subroutine plastic_phenopowerlaw_init !-------------------------------------------------------------------------------------------------- ! optional parameters that need to be defined - prm%c_2 = config%getFloat('twin_b',defaultVal=1.0_pReal) - prm%c_1 = config%getFloat('twin_c',defaultVal=0.0_pReal) - prm%c_4 = config%getFloat('twin_d',defaultVal=0.0_pReal) - prm%c_3 = config%getFloat('twin_e',defaultVal=0.0_pReal) + prm%c_1 = config%getFloat('twin_c',defaultVal=0.0_pReal) + prm%c_2 = config%getFloat('twin_b',defaultVal=1.0_pReal) + prm%c_3 = config%getFloat('twin_e',defaultVal=0.0_pReal) + prm%c_4 = config%getFloat('twin_d',defaultVal=0.0_pReal) prm%aTolResistance = config%getFloat('atol_resistance',defaultVal=1.0_pReal) prm%aTolShear = config%getFloat('atol_shear', defaultVal=1.0e-6_pReal) From 243e884a67c682b3b72c95784689c04ccc752369 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 3 Feb 2020 21:33:40 +0100 Subject: [PATCH 18/20] new test Marc_elementLib should not fail randomly anymore --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index 6f87ed87e..64432754c 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 6f87ed87eddeea4a63827647256fd68328907fca +Subproject commit 64432754ce3c590c882cf4987695539cee524ee8 From 3cadf60d30518cf666d76c4ad3dda186ab2924da Mon Sep 17 00:00:00 2001 From: Test User Date: Wed, 5 Feb 2020 14:49:01 +0100 Subject: [PATCH 19/20] [skip ci] updated version information after successful test of v2.0.3-1609-g243e884a --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 66536aad6..cbd80d7b3 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-1601-gc433e244 +v2.0.3-1609-g243e884a From 177babecb2660d7a8aa9b402abc8a3ff1a2b2ae3 Mon Sep 17 00:00:00 2001 From: Test User Date: Thu, 6 Feb 2020 12:37:13 +0100 Subject: [PATCH 20/20] [skip ci] updated version information after successful test of v2.0.3-1624-g47109b90 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index cbd80d7b3..7e4465e78 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-1609-g243e884a +v2.0.3-1624-g47109b90