From e249168189db9786fe39f0094322a9399a79fb43 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Mon, 22 Feb 2021 18:32:31 +0100 Subject: [PATCH] modifying storing of orientations; needed for nonlocal, and also to remove the use if ip,el at homogenization level. ip, el should be used only for looping eventually. --- src/homogenization.f90 | 6 +++-- src/homogenization_mechanical.f90 | 10 +++---- src/homogenization_thermal.f90 | 11 ++++---- src/material.f90 | 8 +----- src/phase.f90 | 32 ++++++++++++++--------- src/phase_mechanical.f90 | 29 ++++++++++++++------ src/phase_mechanical_plastic_nonlocal.f90 | 9 ++++--- 7 files changed, 62 insertions(+), 43 deletions(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 72ae2231a..0efe2b9af 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -536,13 +536,15 @@ function damage_nonlocal_getDiffusion(ip,el) damage_nonlocal_getDiffusion integer :: & homog, & - grain + grain, & + ce homog = material_homogenizationAt(el) damage_nonlocal_getDiffusion = 0.0_pReal + ce = (el-1)*discretization_nIPs + ip do grain = 1, homogenization_Nconstituents(homog) damage_nonlocal_getDiffusion = damage_nonlocal_getDiffusion + & - crystallite_push33ToRef(grain,ip,el,lattice_D(1:3,1:3,material_phaseAt(grain,el))) + crystallite_push33ToRef(grain,ce,lattice_D(1:3,1:3,material_phaseAt2(grain,ce))) enddo damage_nonlocal_getDiffusion = & diff --git a/src/homogenization_mechanical.f90 b/src/homogenization_mechanical.f90 index cd9eafdb2..74f3e475b 100644 --- a/src/homogenization_mechanical.f90 +++ b/src/homogenization_mechanical.f90 @@ -153,13 +153,13 @@ module subroutine mechanical_homogenize(dt,ip,el) chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) case (HOMOGENIZATION_NONE_ID) chosenHomogenization - homogenization_P(1:3,1:3,ce) = phase_mechanical_getP(1,ip,el) + homogenization_P(1:3,1:3,ce) = phase_mechanical_getP(1,ce) homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = phase_mechanical_dPdF(dt,1,ip,el) case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ip,el) - Ps(:,:,co) = phase_mechanical_getP(co,ip,el) + Ps(:,:,co) = phase_mechanical_getP(co,ce) enddo call mechanical_isostrain_averageStressAndItsTangent(& homogenization_P(1:3,1:3,ce), & @@ -170,7 +170,7 @@ module subroutine mechanical_homogenize(dt,ip,el) case (HOMOGENIZATION_RGC_ID) chosenHomogenization do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ip,el) - Ps(:,:,co) = phase_mechanical_getP(co,ip,el) + Ps(:,:,co) = phase_mechanical_getP(co,ce) enddo call mechanical_RGC_averageStressAndItsTangent(& homogenization_P(1:3,1:3,ce), & @@ -208,8 +208,8 @@ module function mechanical_updateState(subdt,subF,ce,ip,el) result(doneAndHappy) if (homogenization_type(material_homogenizationAt2(ce)) == HOMOGENIZATION_RGC_ID) then do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(subdt,co,ip,el) - Fs(:,:,co) = phase_mechanical_getF(co,ip,el) - Ps(:,:,co) = phase_mechanical_getP(co,ip,el) + Fs(:,:,co) = phase_mechanical_getF(co,ce) + Ps(:,:,co) = phase_mechanical_getP(co,ce) enddo doneAndHappy = mechanical_RGC_updateState(Ps,Fs,subF,subdt,dPdFs,ce) else diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90 index 67027f52a..a02f798cd 100644 --- a/src/homogenization_thermal.f90 +++ b/src/homogenization_thermal.f90 @@ -107,16 +107,17 @@ module function thermal_conduction_getConductivity(ip,el) result(K) real(pReal), dimension(3,3) :: K integer :: & - co - + co, & + ce K = 0.0_pReal - do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - K = K + crystallite_push33ToRef(co,ip,el,lattice_K(:,:,material_phaseAt(co,el))) + ce = (el-1)*discretization_nIPs + ip + do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) + K = K + crystallite_push33ToRef(co,ce,lattice_K(:,:,material_phaseAt2(co,ce))) enddo - K = K / real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) + K = K / real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal) end function thermal_conduction_getConductivity diff --git a/src/material.f90 b/src/material.f90 index 3656f0f0a..c7d24b6f0 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -39,9 +39,6 @@ module material integer, dimension(:,:,:), allocatable, public, protected :: & ! (constituent,IP,elem) material_phaseMemberAt !< position of the element within its phase instance - type(Rotation), dimension(:,:,:), allocatable, public, protected :: & - material_orientation0 !< initial orientation of each grain,IP,element - public :: & material_init @@ -125,8 +122,6 @@ subroutine material_parseMaterial allocate(material_phaseAt2(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0) allocate(material_phaseMemberAt2(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0) - allocate(material_orientation0(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems)) - do el = 1, discretization_Nelems material => materials%get(discretization_materialAt(el)) constituents => material%get('constituents') @@ -135,7 +130,7 @@ subroutine material_parseMaterial do ip = 1, discretization_nIPs ce = (el-1)*discretization_nIPs + ip counterHomogenization(material_homogenizationAt(el)) = counterHomogenization(material_homogenizationAt(el)) + 1 - material_homogenizationMemberAt(ip,el) = counterHomogenization(material_homogenizationAt(el)) + material_homogenizationMemberAt(ip,el) = counterHomogenization(material_homogenizationAt(el)) material_homogenizationAt2(ce) = material_homogenizationAt(el) material_homogenizationMemberAt2(ce) = material_homogenizationMemberAt(ip,el) enddo @@ -153,7 +148,6 @@ subroutine material_parseMaterial material_phaseAt2(co,ce) = material_phaseAt(co,el) material_phaseMemberAt2(co,ce) = material_phaseMemberAt(co,ip,el) - call material_orientation0(co,ip,el)%fromQuaternion(constituent%get_asFloats('O',requiredSize=4)) ! should be done in crystallite enddo enddo diff --git a/src/phase.f90 b/src/phase.f90 index c5de6d658..6bea38161 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -19,6 +19,9 @@ module phase implicit none private + type(Rotation), dimension(:,:,:), allocatable, protected :: & + material_orientation0 !< initial orientation of each grain,IP,element + type(rotation), dimension(:,:,:), allocatable :: & crystallite_orientation !< current orientation @@ -77,8 +80,8 @@ module phase interface ! == cleaned:begin ================================================================================= - module subroutine mechanical_init(phases) - class(tNode), pointer :: phases + module subroutine mechanical_init(materials,phases) + class(tNode), pointer :: materials,phases end subroutine mechanical_init module subroutine damage_init @@ -147,8 +150,8 @@ module phase real(pReal), dimension(3,3) :: L_p end function mechanical_L_p - module function phase_mechanical_getF(co,ip,el) result(F) - integer, intent(in) :: co, ip, el + module function phase_mechanical_getF(co,ce) result(F) + integer, intent(in) :: co, ce real(pReal), dimension(3,3) :: F end function phase_mechanical_getF @@ -157,8 +160,8 @@ module phase real(pReal), dimension(3,3) :: F_e end function mechanical_F_e - module function phase_mechanical_getP(co,ip,el) result(P) - integer, intent(in) :: co, ip, el + module function phase_mechanical_getP(co,ce) result(P) + integer, intent(in) :: co, ce real(pReal), dimension(3,3) :: P end function phase_mechanical_getP @@ -342,6 +345,7 @@ subroutine phase_init so !< counter in source loop class (tNode), pointer :: & debug_constitutive, & + materials, & phases @@ -356,9 +360,10 @@ subroutine phase_init debugConstitutive%grain = config_debug%get_asInt('grain',defaultVal = 1) + materials => config_material%get('material') phases => config_material%get('phase') - call mechanical_init(phases) + call mechanical_init(materials,phases) call damage_init call thermal_init(phases) @@ -624,19 +629,20 @@ end subroutine crystallite_orientations !-------------------------------------------------------------------------------------------------- !> @brief Map 2nd order tensor to reference config !-------------------------------------------------------------------------------------------------- -function crystallite_push33ToRef(co,ip,el, tensor33) +function crystallite_push33ToRef(co,ce, tensor33) real(pReal), dimension(3,3), intent(in) :: tensor33 integer, intent(in):: & - el, & - ip, & - co + co, & + ce real(pReal), dimension(3,3) :: crystallite_push33ToRef real(pReal), dimension(3,3) :: T + integer :: ph, me - - T = matmul(material_orientation0(co,ip,el)%asMatrix(),transpose(math_inv33(phase_mechanical_getF(co,ip,el)))) ! ToDo: initial orientation correct? + ph = material_phaseAt2(co,ce) + me = material_phaseMemberAt2(co,ce) + T = matmul(material_orientation0(co,ph,me)%asMatrix(),transpose(math_inv33(phase_mechanical_getF(co,ce)))) ! ToDo: initial orientation correct? crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T)) diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 24425bf0a..383bf07a8 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -177,21 +177,26 @@ contains !> @brief Initialize mechanical field related constitutive models !> @details Initialize elasticity, plasticity and stiffness degradation models. !-------------------------------------------------------------------------------------------------- -module subroutine mechanical_init(phases) +module subroutine mechanical_init(materials,phases) class(tNode), pointer :: & + materials, & phases integer :: & el, & ip, & co, & + ce, & ph, & me, & stiffDegradationCtr, & Nconstituents class(tNode), pointer :: & num_crystallite, & + material, & + constituents, & + constituent, & phase, & mech, & elastic, & @@ -221,6 +226,8 @@ module subroutine mechanical_init(phases) allocate(phase_mechanical_P(phases%length)) allocate(phase_mechanical_S0(phases%length)) + allocate(material_orientation0(homogenization_maxNconstituents,phases%length,maxval(material_phaseMemberAt))) + do ph = 1, phases%length Nconstituents = count(material_phaseAt == ph) * discretization_nIPs @@ -271,14 +278,20 @@ module subroutine mechanical_init(phases) enddo endif + !$OMP PARALLEL DO PRIVATE(ph,me) do el = 1, size(material_phaseMemberAt,3); do ip = 1, size(material_phaseMemberAt,2) do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) + material => materials%get(discretization_materialAt(el)) + constituents => material%get('constituents') + constituent => constituents%get(co) ph = material_phaseAt(co,el) me = material_phaseMemberAt(co,ip,el) - phase_mechanical_Fp0(ph)%data(1:3,1:3,me) = material_orientation0(co,ip,el)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005) + call material_orientation0(co,ph,me)%fromQuaternion(constituent%get_asFloats('O',requiredSize=4)) + + phase_mechanical_Fp0(ph)%data(1:3,1:3,me) = material_orientation0(co,ph,me)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005) phase_mechanical_Fp0(ph)%data(1:3,1:3,me) = phase_mechanical_Fp0(ph)%data(1:3,1:3,me) & / math_det33(phase_mechanical_Fp0(ph)%data(1:3,1:3,me))**(1.0_pReal/3.0_pReal) phase_mechanical_Fi0(ph)%data(1:3,1:3,me) = math_I3 @@ -1440,13 +1453,13 @@ end function mechanical_L_p !---------------------------------------------------------------------------------------------- !< @brief Get deformation gradient (for use by homogenization) !---------------------------------------------------------------------------------------------- -module function phase_mechanical_getF(co,ip,el) result(F) +module function phase_mechanical_getF(co,ce) result(F) - integer, intent(in) :: co, ip, el + integer, intent(in) :: co, ce real(pReal), dimension(3,3) :: F - F = phase_mechanical_F(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el)) + F = phase_mechanical_F(material_phaseAt2(co,ce))%data(1:3,1:3,material_phaseMemberAt2(co,ce)) end function phase_mechanical_getF @@ -1469,13 +1482,13 @@ end function mechanical_F_e !---------------------------------------------------------------------------------------------- !< @brief Get second Piola-Kichhoff stress (for use by homogenization) !---------------------------------------------------------------------------------------------- -module function phase_mechanical_getP(co,ip,el) result(P) +module function phase_mechanical_getP(co,ce) result(P) - integer, intent(in) :: co, ip, el + integer, intent(in) :: co, ce real(pReal), dimension(3,3) :: P - P = phase_mechanical_P(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el)) + P = phase_mechanical_P(material_phaseAt2(co,ce))%data(1:3,1:3,material_phaseMemberAt2(co,ce)) end function phase_mechanical_getP diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index fbfcfa5af..601618425 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -1403,8 +1403,10 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) integer :: & n, & ! neighbor index + me, & neighbor_e, & ! element index of my neighbor neighbor_i, & ! integration point index of my neighbor + neighbor_me, & neighbor_phase, & ns, & ! number of active slip systems s1, & ! slip system index (me) @@ -1422,6 +1424,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) associate(prm => param(ph)) ns = prm%sum_N_sl + me = material_phaseMemberAt(1,i,e) !*** start out fully compatible my_compatibility = 0.0_pReal forall(s1 = 1:ns) my_compatibility(:,s1,s1,:) = 1.0_pReal @@ -1429,7 +1432,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) neighbors: do n = 1,nIPneighbors neighbor_e = IPneighborhood(1,n,i,e) neighbor_i = IPneighborhood(2,n,i,e) - + neighbor_me = material_phaseMemberAt(1,neighbor_i,neighbor_e) neighbor_phase = material_phaseAt(1,neighbor_e) if (neighbor_e <= 0 .or. neighbor_i <= 0) then @@ -1447,8 +1450,8 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) elseif (prm%chi_GB >= 0.0_pReal) then !* GRAIN BOUNDARY ! !* fixed transmissivity for adjacent ips with different texture (only if explicitly given in material.config) - if (any(dNeq(material_orientation0(1,i,e)%asQuaternion(), & - material_orientation0(1,neighbor_i,neighbor_e)%asQuaternion())) .and. & + if (any(dNeq(material_orientation0(1,ph,me)%asQuaternion(), & + material_orientation0(1,neighbor_phase,neighbor_me)%asQuaternion())) .and. & (.not. phase_localPlasticity(neighbor_phase))) & forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_GB) else