From 4f4adf7d6896584e6e4b39b63ee9bbb8dece0a0d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 27 Jan 2021 07:10:53 +0100 Subject: [PATCH] sorting according to physics --- src/homogenization.f90 | 2 - src/phase.f90 | 84 +------------ src/phase_mechanics.f90 | 141 +++++---------------- src/phase_mechanics_eigendeformation.f90 | 150 +++++++++++++++++++++++ src/phase_mechanics_plastic.f90 | 27 ++++ src/phase_mechanics_plastic_nonlocal.f90 | 1 - src/phase_thermal.f90 | 2 - 7 files changed, 209 insertions(+), 198 deletions(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 4738b3ad8..2255e83b0 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -145,8 +145,6 @@ module homogenization integer, intent(in) :: & ip, & !< integration point number el !< element number - integer :: & - co real(pReal) :: M end function damage_nonlocal_getMobility diff --git a/src/phase.f90 b/src/phase.f90 index d532e4596..c3696eda7 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -15,6 +15,7 @@ module phase use discretization use parallelization use HDF5_utilities + implicit none private @@ -255,8 +256,6 @@ module phase TDot end subroutine constitutive_thermal_getRate - - module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e) integer, intent(in) :: & instance, & @@ -266,54 +265,6 @@ module phase orientation !< crystal orientation end subroutine plastic_nonlocal_updateCompatibility - - module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,me) - 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, & - me - end subroutine plastic_isotropic_LiAndItsTangent - - module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el) - integer, intent(in) :: & - co, & !< grain number - ip, & !< integration point number - el !< element number - real(pReal), intent(in), dimension(3,3) :: & - S - real(pReal), intent(out), dimension(3,3) :: & - Ld !< damage velocity gradient - real(pReal), intent(out), dimension(3,3,3,3) :: & - dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) - end subroutine kinematics_cleavage_opening_LiAndItsTangent - - module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el) - integer, intent(in) :: & - co, & !< grain number - ip, & !< integration point number - el !< element number - real(pReal), intent(in), dimension(3,3) :: & - S - real(pReal), intent(out), dimension(3,3) :: & - Ld !< damage velocity gradient - real(pReal), intent(out), dimension(3,3,3,3) :: & - dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) - end subroutine kinematics_slipplane_opening_LiAndItsTangent - - module subroutine thermalexpansion_LiAndItsTangent(Li, dLi_dTstar, ph,me) - integer, intent(in) :: ph, me - !< element number - real(pReal), intent(out), dimension(3,3) :: & - Li !< thermal velocity gradient - real(pReal), intent(out), dimension(3,3,3,3) :: & - dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero) - end subroutine thermalexpansion_LiAndItsTangent - module subroutine plastic_dependentState(co,ip,el) integer, intent(in) :: & co, & !< component-ID of integration point @@ -350,7 +301,6 @@ module phase constitutive_forward, & constitutive_restore, & plastic_nonlocal_updateCompatibility, & - kinematics_active, & converged, & crystallite_init, & crystallite_stress, & @@ -424,38 +374,6 @@ subroutine constitutive_init end subroutine constitutive_init - -!-------------------------------------------------------------------------------------------------- -!> @brief checks if a kinematic mechanism is active or not -!-------------------------------------------------------------------------------------------------- -function kinematics_active(kinematics_label,kinematics_length) result(active_kinematics) - - character(len=*), intent(in) :: kinematics_label !< name of kinematic mechanism - integer, intent(in) :: kinematics_length !< max. number of kinematics in system - logical, dimension(:,:), allocatable :: active_kinematics - - class(tNode), pointer :: & - phases, & - phase, & - kinematics, & - kinematics_type - integer :: p,k - - phases => config_material%get('phase') - allocate(active_kinematics(kinematics_length,phases%length), source = .false. ) - do p = 1, phases%length - phase => phases%get(p) - kinematics => phase%get('kinematics',defaultVal=emptyList) - do k = 1, kinematics%length - kinematics_type => kinematics%get(k) - if(kinematics_type%get_asString('type') == kinematics_label) active_kinematics(k,p) = .true. - enddo - enddo - - -end function kinematics_active - - !-------------------------------------------------------------------------------------------------- !> @brief Allocate the components of the state structure for a given phase !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_mechanics.f90 b/src/phase_mechanics.f90 index fd8b9331b..dab9ddc88 100644 --- a/src/phase_mechanics.f90 +++ b/src/phase_mechanics.f90 @@ -57,7 +57,17 @@ submodule(phase) mechanics end subroutine plastic_init - + module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,me) + 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, & + me + end subroutine plastic_isotropic_LiAndItsTangent module function plastic_dotState(subdt,co,ip,el,ph,me) result(broken) @@ -85,6 +95,25 @@ submodule(phase) mechanics end function plastic_deltaState + module subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & + S, Fi, co, ip, el) + + integer, intent(in) :: & + co, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + S !< 2nd Piola-Kirchhoff stress + real(pReal), intent(in), dimension(3,3) :: & + Fi !< intermediate deformation gradient + real(pReal), intent(out), dimension(3,3) :: & + Li !< intermediate velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLi_dS, & !< derivative of Li with respect to S + dLi_dFi + + end subroutine constitutive_LiAndItsTangents + module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & S, Fi, co, ip, el) @@ -173,8 +202,7 @@ module subroutine mech_init(phases) phase, & mech, & elastic, & - stiffDegradation, & - kinematics + stiffDegradation print'(/,a)', ' <<<+- constitutive_mech init -+>>>' @@ -318,31 +346,7 @@ module subroutine mech_init(phases) end subroutine mech_init -!-------------------------------------------------------------------------------------------------- -!> @brief checks if a plastic module is active or not -!-------------------------------------------------------------------------------------------------- -function plastic_active(plastic_label) result(active_plastic) - character(len=*), intent(in) :: plastic_label !< type of plasticity model - logical, dimension(:), allocatable :: active_plastic - - class(tNode), pointer :: & - phases, & - phase, & - mech, & - pl - integer :: ph - - phases => config_material%get('phase') - allocate(active_plastic(phases%length), source = .false. ) - do ph = 1, phases%length - phase => phases%get(ph) - mech => phase%get('mechanics') - pl => mech%get('plasticity') - if(pl%get_asString('type') == plastic_label) active_plastic(ph) = .true. - enddo - -end function plastic_active !-------------------------------------------------------------------------------------------------- @@ -486,7 +490,6 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken) ierr, & ! error indicator for LAPACK o, & p, & - m, & ph, & me, & jacoCounterLp, & @@ -1503,86 +1506,4 @@ module subroutine constitutive_mech_setF(F,co,ip,el) end subroutine constitutive_mech_setF -!-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the velocity gradient -! ToDo: MD: S is Mi? -!-------------------------------------------------------------------------------------------------- -subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & - S, Fi, co, ip, el) - - integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element - real(pReal), intent(in), dimension(3,3) :: & - S !< 2nd Piola-Kirchhoff stress - real(pReal), intent(in), dimension(3,3) :: & - Fi !< intermediate deformation gradient - real(pReal), intent(out), dimension(3,3) :: & - Li !< intermediate velocity gradient - real(pReal), intent(out), dimension(3,3,3,3) :: & - dLi_dS, & !< derivative of Li with respect to S - dLi_dFi - - real(pReal), dimension(3,3) :: & - my_Li, & !< intermediate velocity gradient - FiInv, & - temp_33 - real(pReal), dimension(3,3,3,3) :: & - my_dLi_dS - real(pReal) :: & - detFi - integer :: & - k, i, j, & - instance, of, me, ph - - Li = 0.0_pReal - dLi_dS = 0.0_pReal - dLi_dFi = 0.0_pReal - - plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) - case (PLASTICITY_isotropic_ID) plasticityType - of = material_phasememberAt(co,ip,el) - instance = phase_plasticityInstance(material_phaseAt(co,el)) - call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,instance,of) - case default plasticityType - my_Li = 0.0_pReal - my_dLi_dS = 0.0_pReal - end select plasticityType - - Li = Li + my_Li - dLi_dS = dLi_dS + my_dLi_dS - - KinematicsLoop: do k = 1, phase_Nkinematics(material_phaseAt(co,el)) - kinematicsType: select case (phase_kinematics(k,material_phaseAt(co,el))) - case (KINEMATICS_cleavage_opening_ID) kinematicsType - call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el) - case (KINEMATICS_slipplane_opening_ID) kinematicsType - call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el) - case (KINEMATICS_thermal_expansion_ID) kinematicsType - me = material_phaseMemberAt(co,ip,el) - ph = material_phaseAt(co,el) - call thermalexpansion_LiAndItsTangent(my_Li, my_dLi_dS, ph,me) - case default kinematicsType - my_Li = 0.0_pReal - my_dLi_dS = 0.0_pReal - end select kinematicsType - Li = Li + my_Li - dLi_dS = dLi_dS + my_dLi_dS - enddo KinematicsLoop - - FiInv = math_inv33(Fi) - detFi = math_det33(Fi) - Li = matmul(matmul(Fi,Li),FiInv)*detFi !< push forward to intermediate configuration - temp_33 = matmul(FiInv,Li) - - do i = 1,3; do j = 1,3 - dLi_dS(1:3,1:3,i,j) = matmul(matmul(Fi,dLi_dS(1:3,1:3,i,j)),FiInv)*detFi - dLi_dFi(1:3,1:3,i,j) = dLi_dFi(1:3,1:3,i,j) + Li*FiInv(j,i) - dLi_dFi(1:3,i,1:3,j) = dLi_dFi(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i) - enddo; enddo - -end subroutine constitutive_LiAndItsTangents - end submodule mechanics - diff --git a/src/phase_mechanics_eigendeformation.f90 b/src/phase_mechanics_eigendeformation.f90 index 98bc8481f..d1f581d5a 100644 --- a/src/phase_mechanics_eigendeformation.f90 +++ b/src/phase_mechanics_eigendeformation.f90 @@ -15,6 +15,42 @@ submodule(phase:mechanics) eigendeformation integer, intent(in) :: kinematics_length logical, dimension(:,:), allocatable :: myKinematics end function kinematics_thermal_expansion_init + + module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el) + integer, intent(in) :: & + co, & !< grain number + ip, & !< integration point number + el !< element number + real(pReal), intent(in), dimension(3,3) :: & + S + real(pReal), intent(out), dimension(3,3) :: & + Ld !< damage velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) + end subroutine kinematics_cleavage_opening_LiAndItsTangent + + module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el) + integer, intent(in) :: & + co, & !< grain number + ip, & !< integration point number + el !< element number + real(pReal), intent(in), dimension(3,3) :: & + S + real(pReal), intent(out), dimension(3,3) :: & + Ld !< damage velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) + end subroutine kinematics_slipplane_opening_LiAndItsTangent + + module subroutine thermalexpansion_LiAndItsTangent(Li, dLi_dTstar, ph,me) + integer, intent(in) :: ph, me + !< element number + real(pReal), intent(out), dimension(3,3) :: & + Li !< thermal velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero) + end subroutine thermalexpansion_LiAndItsTangent + end interface @@ -54,4 +90,118 @@ module subroutine eigendeformation_init(phases) end subroutine eigendeformation_init + +!-------------------------------------------------------------------------------------------------- +!> @brief checks if a kinematic mechanism is active or not +!-------------------------------------------------------------------------------------------------- +function kinematics_active(kinematics_label,kinematics_length) result(active_kinematics) + + character(len=*), intent(in) :: kinematics_label !< name of kinematic mechanism + integer, intent(in) :: kinematics_length !< max. number of kinematics in system + logical, dimension(:,:), allocatable :: active_kinematics + + class(tNode), pointer :: & + phases, & + phase, & + kinematics, & + kinematics_type + integer :: p,k + + phases => config_material%get('phase') + allocate(active_kinematics(kinematics_length,phases%length), source = .false. ) + do p = 1, phases%length + phase => phases%get(p) + kinematics => phase%get('kinematics',defaultVal=emptyList) + do k = 1, kinematics%length + kinematics_type => kinematics%get(k) + if(kinematics_type%get_asString('type') == kinematics_label) active_kinematics(k,p) = .true. + enddo + enddo + + +end function kinematics_active + + +!-------------------------------------------------------------------------------------------------- +!> @brief contains the constitutive equation for calculating the velocity gradient +! ToDo: MD: S is Mi? +!-------------------------------------------------------------------------------------------------- +module subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & + S, Fi, co, ip, el) + + integer, intent(in) :: & + co, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + S !< 2nd Piola-Kirchhoff stress + real(pReal), intent(in), dimension(3,3) :: & + Fi !< intermediate deformation gradient + real(pReal), intent(out), dimension(3,3) :: & + Li !< intermediate velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLi_dS, & !< derivative of Li with respect to S + dLi_dFi + + real(pReal), dimension(3,3) :: & + my_Li, & !< intermediate velocity gradient + FiInv, & + temp_33 + real(pReal), dimension(3,3,3,3) :: & + my_dLi_dS + real(pReal) :: & + detFi + integer :: & + k, i, j, & + instance, of, me, ph + + Li = 0.0_pReal + dLi_dS = 0.0_pReal + dLi_dFi = 0.0_pReal + + plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) + case (PLASTICITY_isotropic_ID) plasticityType + of = material_phasememberAt(co,ip,el) + instance = phase_plasticityInstance(material_phaseAt(co,el)) + call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,instance,of) + case default plasticityType + my_Li = 0.0_pReal + my_dLi_dS = 0.0_pReal + end select plasticityType + + Li = Li + my_Li + dLi_dS = dLi_dS + my_dLi_dS + + KinematicsLoop: do k = 1, phase_Nkinematics(material_phaseAt(co,el)) + kinematicsType: select case (phase_kinematics(k,material_phaseAt(co,el))) + case (KINEMATICS_cleavage_opening_ID) kinematicsType + call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el) + case (KINEMATICS_slipplane_opening_ID) kinematicsType + call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el) + case (KINEMATICS_thermal_expansion_ID) kinematicsType + me = material_phaseMemberAt(co,ip,el) + ph = material_phaseAt(co,el) + call thermalexpansion_LiAndItsTangent(my_Li, my_dLi_dS, ph,me) + case default kinematicsType + my_Li = 0.0_pReal + my_dLi_dS = 0.0_pReal + end select kinematicsType + Li = Li + my_Li + dLi_dS = dLi_dS + my_dLi_dS + enddo KinematicsLoop + + FiInv = math_inv33(Fi) + detFi = math_det33(Fi) + Li = matmul(matmul(Fi,Li),FiInv)*detFi !< push forward to intermediate configuration + temp_33 = matmul(FiInv,Li) + + do i = 1,3; do j = 1,3 + dLi_dS(1:3,1:3,i,j) = matmul(matmul(Fi,dLi_dS(1:3,1:3,i,j)),FiInv)*detFi + dLi_dFi(1:3,1:3,i,j) = dLi_dFi(1:3,1:3,i,j) + Li*FiInv(j,i) + dLi_dFi(1:3,i,1:3,j) = dLi_dFi(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i) + enddo; enddo + +end subroutine constitutive_LiAndItsTangents + + end submodule eigendeformation diff --git a/src/phase_mechanics_plastic.f90 b/src/phase_mechanics_plastic.f90 index 9b8334c99..b03683061 100644 --- a/src/phase_mechanics_plastic.f90 +++ b/src/phase_mechanics_plastic.f90 @@ -439,4 +439,31 @@ module function plastic_deltaState(co, ip, el, ph, of) result(broken) end function plastic_deltaState + +!-------------------------------------------------------------------------------------------------- +!> @brief checks if a plastic module is active or not +!-------------------------------------------------------------------------------------------------- +function plastic_active(plastic_label) result(active_plastic) + + character(len=*), intent(in) :: plastic_label !< type of plasticity model + logical, dimension(:), allocatable :: active_plastic + + class(tNode), pointer :: & + phases, & + phase, & + mech, & + pl + integer :: ph + + phases => config_material%get('phase') + allocate(active_plastic(phases%length), source = .false. ) + do ph = 1, phases%length + phase => phases%get(ph) + mech => phase%get('mechanics') + pl => mech%get('plasticity') + if(pl%get_asString('type') == plastic_label) active_plastic(ph) = .true. + enddo + +end function plastic_active + end submodule plastic diff --git a/src/phase_mechanics_plastic_nonlocal.f90 b/src/phase_mechanics_plastic_nonlocal.f90 index d0ed01cf0..2f566c515 100644 --- a/src/phase_mechanics_plastic_nonlocal.f90 +++ b/src/phase_mechanics_plastic_nonlocal.f90 @@ -988,7 +988,6 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & el !< current element number integer :: & - ph, & ns, & !< short notation for the total number of active slip systems c, & !< character of dislocation t, & !< type of dislocation diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90 index 621bb4f4a..04c35a311 100644 --- a/src/phase_thermal.f90 +++ b/src/phase_thermal.f90 @@ -187,8 +187,6 @@ module function thermal_stress(Delta_t,ph,me) result(converged_) integer, intent(in) :: ph, me logical :: converged_ - integer :: so - converged_ = .not. integrateThermalState(Delta_t,ph,me)