From 912a21f5b6fe1098d1ca5be7eeb4ed1072ed0b9a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 26 Jan 2021 01:11:32 +0100 Subject: [PATCH] modularizing --- src/phase_mechanics.f90 | 169 ++---------------- src/phase_mechanics_plastic.f90 | 162 +++++++++++++++++ src/phase_mechanics_plastic_dislotungsten.f90 | 4 +- src/phase_mechanics_plastic_dislotwin.f90 | 4 +- src/phase_mechanics_plastic_isotropic.f90 | 4 +- src/phase_mechanics_plastic_kinehardening.f90 | 4 +- src/phase_mechanics_plastic_none.f90 | 4 +- src/phase_mechanics_plastic_nonlocal.f90 | 4 +- src/phase_mechanics_plastic_phenopowerlaw.f90 | 4 +- 9 files changed, 192 insertions(+), 167 deletions(-) create mode 100644 src/phase_mechanics_plastic.f90 diff --git a/src/phase_mechanics.f90 b/src/phase_mechanics.f90 index 53986191b..fe0a95fbb 100644 --- a/src/phase_mechanics.f90 +++ b/src/phase_mechanics.f90 @@ -83,93 +83,6 @@ submodule(constitutive) constitutive_mech myPlasticity end function plastic_nonlocal_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 - - pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) - real(pReal), dimension(3,3), intent(out) :: & - Lp !< plastic velocity gradient - real(pReal), dimension(3,3,3,3), intent(out) :: & - dLp_dMp !< derivative of Lp with respect to the Mandel stress - - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - end subroutine plastic_kinehardening_LpAndItsTangent - - module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) - real(pReal), dimension(3,3), intent(out) :: & - Lp !< plastic velocity gradient - real(pReal), dimension(3,3,3,3), intent(out) :: & - dLp_dMp !< derivative of Lp with respect to the Mandel stress - - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - real(pReal), intent(in) :: & - T - integer, intent(in) :: & - instance, & - of - end subroutine plastic_dislotwin_LpAndItsTangent - - pure module subroutine plastic_dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) - real(pReal), dimension(3,3), intent(out) :: & - Lp !< plastic velocity gradient - real(pReal), dimension(3,3,3,3), intent(out) :: & - dLp_dMp !< derivative of Lp with respect to the Mandel stress - - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - real(pReal), intent(in) :: & - T - integer, intent(in) :: & - instance, & - of - end subroutine plastic_dislotungsten_LpAndItsTangent - - module subroutine plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp, & - Mp,Temperature,instance,of,ip,el) - real(pReal), dimension(3,3), intent(out) :: & - Lp !< plastic velocity gradient - real(pReal), dimension(3,3,3,3), intent(out) :: & - dLp_dMp !< derivative of Lp with respect to the Mandel stress - - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - real(pReal), intent(in) :: & - Temperature - integer, intent(in) :: & - instance, & - of, & - ip, & !< current integration point - el !< current element number - end subroutine plastic_nonlocal_LpAndItsTangent - module subroutine plastic_isotropic_dotState(Mp,instance,of) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -267,6 +180,22 @@ submodule(constitutive) constitutive_mech ip, & el end subroutine plastic_nonlocal_deltaState +module subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & + S, Fi, co, ip, el) + integer, intent(in) :: & + co, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + S, & !< 2nd Piola-Kirchhoff stress + Fi !< intermediate deformation gradient + real(pReal), intent(out), dimension(3,3) :: & + Lp !< plastic velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLp_dS, & + dLp_dFi !< derivative of Lp with respect to Fi + +end subroutine constitutive_plastic_LpAndItsTangents module function kinematics_cleavage_opening_init(kinematics_length) result(myKinematics) integer, intent(in) :: kinematics_length @@ -621,72 +550,6 @@ module subroutine constitutive_plastic_dependentState(co, ip, el) end subroutine constitutive_plastic_dependentState -!-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the velocity gradient -! ToDo: Discuss whether it makes sense if crystallite handles the configuration conversion, i.e. -! Mp in, dLp_dMp out -!-------------------------------------------------------------------------------------------------- -subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & - S, Fi, co, ip, el) - integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element - real(pReal), intent(in), dimension(3,3) :: & - S, & !< 2nd Piola-Kirchhoff stress - Fi !< intermediate deformation gradient - real(pReal), intent(out), dimension(3,3) :: & - Lp !< plastic velocity gradient - real(pReal), intent(out), dimension(3,3,3,3) :: & - dLp_dS, & - dLp_dFi !< derivative of Lp with respect to Fi - - real(pReal), dimension(3,3,3,3) :: & - dLp_dMp !< derivative of Lp with respect to Mandel stress - real(pReal), dimension(3,3) :: & - Mp !< Mandel stress work conjugate with Lp - integer :: & - i, j, instance, me, ph - - - Mp = matmul(matmul(transpose(Fi),Fi),S) - me = material_phasememberAt(co,ip,el) - ph = material_phaseAt(co,el) - instance = phase_plasticityInstance(ph) - - plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) - - case (PLASTICITY_NONE_ID) plasticityType - Lp = 0.0_pReal - dLp_dMp = 0.0_pReal - - case (PLASTICITY_ISOTROPIC_ID) plasticityType - call plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me) - - case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType - call plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me) - - case (PLASTICITY_KINEHARDENING_ID) plasticityType - call plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me) - - case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),instance,me,ip,el) - - case (PLASTICITY_DISLOTWIN_ID) plasticityType - call plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),instance,me) - - case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType - call plastic_dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),instance,me) - - end select plasticityType - - do i=1,3; do j=1,3 - dLp_dFi(i,j,1:3,1:3) = matmul(matmul(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + & - matmul(matmul(Fi,dLp_dMp(i,j,1:3,1:3)),S) - dLp_dS(i,j,1:3,1:3) = matmul(matmul(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi) - enddo; enddo - -end subroutine constitutive_plastic_LpAndItsTangents !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_mechanics_plastic.f90 b/src/phase_mechanics_plastic.f90 new file mode 100644 index 000000000..fdcc22027 --- /dev/null +++ b/src/phase_mechanics_plastic.f90 @@ -0,0 +1,162 @@ +submodule(constitutive:constitutive_mech) plastic + + interface + + 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 + + pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) + real(pReal), dimension(3,3), intent(out) :: & + Lp !< plastic velocity gradient + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp !< derivative of Lp with respect to the Mandel stress + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + end subroutine plastic_kinehardening_LpAndItsTangent + + module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) + real(pReal), dimension(3,3), intent(out) :: & + Lp !< plastic velocity gradient + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp !< derivative of Lp with respect to the Mandel stress + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + T + integer, intent(in) :: & + instance, & + of + end subroutine plastic_dislotwin_LpAndItsTangent + + pure module subroutine plastic_dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) + real(pReal), dimension(3,3), intent(out) :: & + Lp !< plastic velocity gradient + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp !< derivative of Lp with respect to the Mandel stress + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + T + integer, intent(in) :: & + instance, & + of + end subroutine plastic_dislotungsten_LpAndItsTangent + + module subroutine plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp, & + Mp,Temperature,instance,of,ip,el) + real(pReal), dimension(3,3), intent(out) :: & + Lp !< plastic velocity gradient + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp !< derivative of Lp with respect to the Mandel stress + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + Temperature + integer, intent(in) :: & + instance, & + of, & + ip, & !< current integration point + el !< current element number + end subroutine plastic_nonlocal_LpAndItsTangent + + end interface + +contains + +!-------------------------------------------------------------------------------------------------- +!> @brief contains the constitutive equation for calculating the velocity gradient +! ToDo: Discuss whether it makes sense if crystallite handles the configuration conversion, i.e. +! Mp in, dLp_dMp out +!-------------------------------------------------------------------------------------------------- +module subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & + S, Fi, co, ip, el) + integer, intent(in) :: & + co, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + S, & !< 2nd Piola-Kirchhoff stress + Fi !< intermediate deformation gradient + real(pReal), intent(out), dimension(3,3) :: & + Lp !< plastic velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLp_dS, & + dLp_dFi !< derivative of Lp with respect to Fi + + real(pReal), dimension(3,3,3,3) :: & + dLp_dMp !< derivative of Lp with respect to Mandel stress + real(pReal), dimension(3,3) :: & + Mp !< Mandel stress work conjugate with Lp + integer :: & + i, j, instance, me, ph + + + Mp = matmul(matmul(transpose(Fi),Fi),S) + me = material_phasememberAt(co,ip,el) + ph = material_phaseAt(co,el) + instance = phase_plasticityInstance(ph) + + plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) + + case (PLASTICITY_NONE_ID) plasticityType + Lp = 0.0_pReal + dLp_dMp = 0.0_pReal + + case (PLASTICITY_ISOTROPIC_ID) plasticityType + call plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me) + + case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType + call plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me) + + case (PLASTICITY_KINEHARDENING_ID) plasticityType + call plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me) + + case (PLASTICITY_NONLOCAL_ID) plasticityType + call plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),instance,me,ip,el) + + case (PLASTICITY_DISLOTWIN_ID) plasticityType + call plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),instance,me) + + case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType + call plastic_dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),instance,me) + + end select plasticityType + + do i=1,3; do j=1,3 + dLp_dFi(i,j,1:3,1:3) = matmul(matmul(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + & + matmul(matmul(Fi,dLp_dMp(i,j,1:3,1:3)),S) + dLp_dS(i,j,1:3,1:3) = matmul(matmul(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi) + enddo; enddo + +end subroutine constitutive_plastic_LpAndItsTangents + +end submodule plastic diff --git a/src/phase_mechanics_plastic_dislotungsten.f90 b/src/phase_mechanics_plastic_dislotungsten.f90 index c39ae5c2b..8b258ab2a 100644 --- a/src/phase_mechanics_plastic_dislotungsten.f90 +++ b/src/phase_mechanics_plastic_dislotungsten.f90 @@ -5,7 +5,7 @@ !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief crystal plasticity model for bcc metals, especially Tungsten !-------------------------------------------------------------------------------------------------- -submodule(constitutive:constitutive_mech) plastic_dislotungsten +submodule(constitutive:plastic) dislotungsten real(pReal), parameter :: & kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin @@ -547,4 +547,4 @@ pure subroutine kinetics(Mp,T,instance,of, & end subroutine kinetics -end submodule plastic_dislotungsten +end submodule dislotungsten diff --git a/src/phase_mechanics_plastic_dislotwin.f90 b/src/phase_mechanics_plastic_dislotwin.f90 index dea84f751..c2b19013e 100644 --- a/src/phase_mechanics_plastic_dislotwin.f90 +++ b/src/phase_mechanics_plastic_dislotwin.f90 @@ -7,7 +7,7 @@ !> @brief material subroutine incoprorating dislocation and twinning physics !> @details to be done !-------------------------------------------------------------------------------------------------- -submodule(constitutive:constitutive_mech) plastic_dislotwin +submodule(constitutive:plastic) dislotwin real(pReal), parameter :: & kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin @@ -1088,4 +1088,4 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,of,& end subroutine kinetics_trans -end submodule plastic_dislotwin +end submodule dislotwin diff --git a/src/phase_mechanics_plastic_isotropic.f90 b/src/phase_mechanics_plastic_isotropic.f90 index b7c5f67c1..eab9b67cc 100644 --- a/src/phase_mechanics_plastic_isotropic.f90 +++ b/src/phase_mechanics_plastic_isotropic.f90 @@ -7,7 +7,7 @@ !! resolving the stress on the slip systems. Will give the response of phenopowerlaw for an !! untextured polycrystal !-------------------------------------------------------------------------------------------------- -submodule(constitutive:constitutive_mech) plastic_isotropic +submodule(constitutive:plastic) isotropic type :: tParameters real(pReal) :: & @@ -348,4 +348,4 @@ module subroutine plastic_isotropic_results(instance,group) end subroutine plastic_isotropic_results -end submodule plastic_isotropic +end submodule isotropic diff --git a/src/phase_mechanics_plastic_kinehardening.f90 b/src/phase_mechanics_plastic_kinehardening.f90 index 8454b28f8..4a92f6567 100644 --- a/src/phase_mechanics_plastic_kinehardening.f90 +++ b/src/phase_mechanics_plastic_kinehardening.f90 @@ -5,7 +5,7 @@ !> @brief Phenomenological crystal plasticity using a power law formulation for the shear rates !! and a Voce-type kinematic hardening rule !-------------------------------------------------------------------------------------------------- -submodule(constitutive:constitutive_mech) plastic_kinehardening +submodule(constitutive:plastic) kinehardening type :: tParameters real(pReal) :: & @@ -474,4 +474,4 @@ pure subroutine kinetics(Mp,instance,of, & end subroutine kinetics -end submodule plastic_kinehardening +end submodule kinehardening diff --git a/src/phase_mechanics_plastic_none.f90 b/src/phase_mechanics_plastic_none.f90 index 27a01fb93..88fa6a5c2 100644 --- a/src/phase_mechanics_plastic_none.f90 +++ b/src/phase_mechanics_plastic_none.f90 @@ -4,7 +4,7 @@ !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief Dummy plasticity for purely elastic material !-------------------------------------------------------------------------------------------------- -submodule(constitutive:constitutive_mech) plastic_none +submodule(constitutive:plastic) none contains @@ -50,4 +50,4 @@ module function plastic_none_init() result(myPlasticity) end function plastic_none_init -end submodule plastic_none +end submodule none diff --git a/src/phase_mechanics_plastic_nonlocal.f90 b/src/phase_mechanics_plastic_nonlocal.f90 index 2244eb7ad..648a6cbe4 100644 --- a/src/phase_mechanics_plastic_nonlocal.f90 +++ b/src/phase_mechanics_plastic_nonlocal.f90 @@ -4,7 +4,7 @@ !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @brief material subroutine for plasticity including dislocation flux !-------------------------------------------------------------------------------------------------- -submodule(constitutive:constitutive_mech) plastic_nonlocal +submodule(constitutive:plastic) nonlocal use geometry_plastic_nonlocal, only: & nIPneighbors => geometry_plastic_nonlocal_nIPneighbors, & IPneighborhood => geometry_plastic_nonlocal_IPneighborhood, & @@ -1834,4 +1834,4 @@ pure function getRho0(instance,of,ip,el) end function getRho0 -end submodule plastic_nonlocal +end submodule nonlocal diff --git a/src/phase_mechanics_plastic_phenopowerlaw.f90 b/src/phase_mechanics_plastic_phenopowerlaw.f90 index 678acad27..39ceb1000 100644 --- a/src/phase_mechanics_plastic_phenopowerlaw.f90 +++ b/src/phase_mechanics_plastic_phenopowerlaw.f90 @@ -4,7 +4,7 @@ !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief phenomenological crystal plasticity formulation using a powerlaw fitting !-------------------------------------------------------------------------------------------------- -submodule(constitutive:constitutive_mech) plastic_phenopowerlaw +submodule(constitutive:plastic) phenopowerlaw type :: tParameters real(pReal) :: & @@ -543,4 +543,4 @@ pure subroutine kinetics_twin(Mp,instance,of,& end subroutine kinetics_twin -end submodule plastic_phenopowerlaw +end submodule phenopowerlaw