From c70e535da83d876d52b636149503a5d5f3365871 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 4 Feb 2022 00:11:08 +0100 Subject: [PATCH 01/13] avoid access by el/ip --- src/phase_mechanical_plastic_nonlocal.f90 | 36 ++++++++++++----------- 1 file changed, 19 insertions(+), 17 deletions(-) diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index a8f56b5c5..81ce4d08e 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -123,8 +123,10 @@ submodule(phase:plastic) nonlocal type :: tNonlocalDependentState real(pReal), allocatable, dimension(:,:) :: & - tau_pass, & - tau_Back + tau_pass, & + tau_Back + real(pReal), allocatable, dimension(:,:,:,:,:) :: & + compatibility end type tNonlocalDependentState type :: tNonlocalState @@ -160,7 +162,6 @@ submodule(phase:plastic) nonlocal state0 type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters - type(tNonlocalDependentState), dimension(:), allocatable :: dependentState contains @@ -489,6 +490,7 @@ module function plastic_nonlocal_init() result(myPlasticity) allocate(dst%tau_pass(prm%sum_N_sl,Nmembers),source=0.0_pReal) allocate(dst%tau_back(prm%sum_N_sl,Nmembers),source=0.0_pReal) + allocate(dst%compatibility(2,maxval(param%sum_N_sl),maxval(param%sum_N_sl),nIPneighbors,Nmembers),source=0.0_pReal) end associate if (Nmembers > 0) call stateInit(ini,ph,Nmembers) @@ -669,7 +671,7 @@ module subroutine nonlocal_dependentState(ph, en, ip, el) - discretization_IPcoords(1:3,el+neighbor_ip-1)) normal_latticeConf = matmul(transpose(invFp), IPareaNormal(1:3,n,ip,el)) if (math_inner(normal_latticeConf,connection_latticeConf(1:3,n)) < 0.0_pReal) & ! neighboring connection points in opposite direction to face normal: must be periodic image - connection_latticeConf(1:3,n) = normal_latticeConf * IPvolume(ip,el)/IParea(n,ip,el) ! instead take the surface normal scaled with the diameter of the cell + connection_latticeConf(1:3,n) = normal_latticeConf * geom(ph)%V_0(en)/IParea(n,ip,el) ! instead take the surface normal scaled with the diameter of the cell else ! local neighbor or different lattice structure or different constitution instance -> use central values instead connection_latticeConf(1:3,n) = 0.0_pReal @@ -1110,7 +1112,7 @@ module subroutine nonlocal_dotState(Mp,timestep, & .or. any(rho(:,dip) + rhoDot(:,9:10) * timestep < -prm%atol_rho)) then #ifdef DEBUG if (debugConstitutive%extensive) then - print'(a,i5,a,i2)', '<< CONST >> evolution rate leads to negative density at el ',el,' ip ',ip + print'(a,i5,a,i2)', '<< CONST >> evolution rate leads to negative density at ph ',ph,' en ',en print'(a)', '<< CONST >> enforcing cutback !!!' end if #endif @@ -1215,14 +1217,14 @@ function rhoDotFlux(timestep,ph,en,ip,el) !*** check CFL (Courant-Friedrichs-Lewy) condition for flux if (any( abs(dot_gamma) > 0.0_pReal & ! any active slip system ... .and. prm%C_CFL * abs(v0) * timestep & - > IPvolume(ip,el) / maxval(IParea(:,ip,el)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here) + > geom(ph)%V_0(en)/ maxval(IParea(:,ip,el)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here) #ifdef DEBUG if (debugConstitutive%extensive) then - print'(a,i5,a,i2)', '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip + print'(a,i5,a,i2)', '<< CONST >> CFL condition not fullfilled at ph ',ph,' en ',en print'(a,e10.3,a,e10.3)', '<< CONST >> velocity is at ', & maxval(abs(v0), abs(dot_gamma) > 0.0_pReal & .and. prm%C_CFL * abs(v0) * timestep & - > IPvolume(ip,el) / maxval(IParea(:,ip,el))), & + > geom(ph)%V_0(en) / maxval(IParea(:,ip,el))), & ' at a timestep of ',timestep print*, '<< CONST >> enforcing cutback !!!' end if @@ -1277,7 +1279,7 @@ function rhoDotFlux(timestep,ph,en,ip,el) !* compatibility if (neighbor_n > 0) then if (phase_plasticity(material_phaseAt(1,neighbor_el)) == PLASTIC_NONLOCAL_ID .and. & - any(compatibility(:,:,:,n,ip,el) > 0.0_pReal)) then + any(dependentState(ph)%compatibility(:,:,:,n,en) > 0.0_pReal)) then forall (s = 1:ns, t = 1:4) neighbor_v0(s,t) = plasticState(np)%state0(iV (s,t,neighbor_ph),no) @@ -1301,12 +1303,12 @@ function rhoDotFlux(timestep,ph,en,ip,el) .and. v0(s,t) * neighbor_v0(s,t) >= 0.0_pReal ) then ! ... only if no sign change in flux density lineLength = neighbor_rhoSgl0(s,t) * neighbor_v0(s,t) & * math_inner(m(1:3,s,t), normal_neighbor2me) * area ! positive line length that wants to enter through this interface - where (compatibility(c,:,s,n,ip,el) > 0.0_pReal) & + where (dependentState(ph)%compatibility(c,:,s,n,en) > 0.0_pReal) & rhoDotFlux(:,t) = rhoDotFlux(1:ns,t) & - + lineLength/IPvolume(ip,el)*compatibility(c,:,s,n,ip,el)**2 ! transferring to equally signed mobile dislocation type - where (compatibility(c,:,s,n,ip,el) < 0.0_pReal) & + + lineLength/geom(ph)%V_0(en)*dependentState(ph)%compatibility(c,:,s,n,en)**2 ! transferring to equally signed mobile dislocation type + where (dependentState(ph)%compatibility(c,:,s,n,en) < 0.0_pReal) & rhoDotFlux(:,topp) = rhoDotFlux(:,topp) & - + lineLength/IPvolume(ip,el)*compatibility(c,:,s,n,ip,el)**2 ! transferring to opposite signed mobile dislocation type + + lineLength/geom(ph)%V_0(en)*dependentState(ph)%compatibility(c,:,s,n,en)**2 ! transferring to opposite signed mobile dislocation type end if end do @@ -1335,15 +1337,15 @@ function rhoDotFlux(timestep,ph,en,ip,el) c = (t + 1) / 2 if (v0(s,t) * math_inner(m(1:3,s,t), normal_me2neighbor) > 0.0_pReal ) then ! flux from en to my neighbor == leaving flux for en (might also be a pure flux from my mobile density to dead density if interface not at all transmissive) if (v0(s,t) * neighbor_v0(s,t) >= 0.0_pReal) then ! no sign change in flux density - transmissivity = sum(compatibility(c,:,s,n,ip,el)**2) ! overall transmissivity from this slip system to my neighbor + transmissivity = sum(dependentState(ph)%compatibility(c,:,s,n,en)**2) ! overall transmissivity from this slip system to my neighbor else ! sign change in flux density means sign change in stress which does not allow for dislocations to arive at the neighbor transmissivity = 0.0_pReal end if lineLength = my_rhoSgl0(s,t) * v0(s,t) & * math_inner(m(1:3,s,t), normal_me2neighbor) * area ! positive line length of mobiles that wants to leave through this interface - rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / IPvolume(ip,el) ! subtract dislocation flux from current type + rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / geom(ph)%V_0(en) ! subtract dislocation flux from current type rhoDotFlux(s,t+4) = rhoDotFlux(s,t+4) & - + lineLength / IPvolume(ip,el) * (1.0_pReal - transmissivity) & + + lineLength / geom(ph)%V_0(en) * (1.0_pReal - transmissivity) & * sign(1.0_pReal, v0(s,t)) ! dislocation flux that is not able to leave through interface (because of low transmissivity) will remain as immobile single density at the material point end if end do @@ -1465,7 +1467,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) end do neighbors - compatibility(:,:,:,:,i,e) = my_compatibility + dependentState(ph)%compatibility(:,:,:,:,material_phaseEntry(1,(e-1)*discretization_nIPs + i)) = my_compatibility end associate From 8b12814eec5678221d7b0e2535a9d5f53a105266 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 4 Feb 2022 06:42:14 +0100 Subject: [PATCH 02/13] common variable names --- src/phase.f90 | 6 +++--- src/phase_mechanical_plastic_nonlocal.f90 | 14 +++++++------- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/phase.f90 b/src/phase.f90 index b7bace81a..82bc44bfd 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -264,11 +264,11 @@ module phase real(pReal) :: f end function phase_f_T - module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) + module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,ip,el) integer, intent(in) :: & ph, & - i, & - e + ip, & + el type(tRotationContainer), dimension(:), intent(in) :: orientation end subroutine plastic_nonlocal_updateCompatibility diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 81ce4d08e..bb36aee06 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -1366,14 +1366,14 @@ end function rhoDotFlux ! 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. !-------------------------------------------------------------------------------------------------- -module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) +module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,ip,el) type(tRotationContainer), dimension(:), intent(in) :: & orientation ! crystal orientation integer, intent(in) :: & ph, & - i, & - e + ip, & + el integer :: & n, & ! neighbor index @@ -1399,14 +1399,14 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) associate(prm => param(ph)) ns = prm%sum_N_sl - en = material_phaseMemberAt(1,i,e) + en = material_phaseMemberAt(1,ip,el) !*** start out fully compatible my_compatibility = 0.0_pReal forall(s1 = 1:ns) my_compatibility(:,s1,s1,:) = 1.0_pReal neighbors: do n = 1,nIPneighbors - neighbor_e = IPneighborhood(1,n,i,e) - neighbor_i = IPneighborhood(2,n,i,e) + neighbor_e = IPneighborhood(1,n,ip,el) + neighbor_i = IPneighborhood(2,n,ip,el) neighbor_me = material_phaseMemberAt(1,neighbor_i,neighbor_e) neighbor_phase = material_phaseAt(1,neighbor_e) @@ -1467,7 +1467,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) end do neighbors - dependentState(ph)%compatibility(:,:,:,:,material_phaseEntry(1,(e-1)*discretization_nIPs + i)) = my_compatibility + dependentState(ph)%compatibility(:,:,:,:,material_phaseEntry(1,(el-1)*discretization_nIPs + ip)) = my_compatibility end associate From 900ef0a2c96544d5d72246e86f0f00e3ee6ad93f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 4 Feb 2022 07:33:12 +0100 Subject: [PATCH 03/13] store IP neighborhood for ph/en access --- src/phase_mechanical_plastic_nonlocal.f90 | 63 ++++++++++++++--------- 1 file changed, 40 insertions(+), 23 deletions(-) diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index bb36aee06..78271b92f 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -15,6 +15,9 @@ submodule(phase:plastic) nonlocal type :: tGeometry real(pReal), dimension(:), allocatable :: V_0 + integer, dimension(:,:,:), allocatable :: IPneighborhood + real(pReal), dimension(:,:), allocatable :: IParea + real(pReal), dimension(:,:,:), allocatable :: IPareaNormal end type tGeometry type(tGeometry), dimension(:), allocatable :: geom @@ -407,6 +410,9 @@ module function plastic_nonlocal_init() result(myPlasticity) call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState,0) ! ToDo: state structure does not follow convention allocate(geom(ph)%V_0(Nmembers)) + allocate(geom(ph)%IPneighborhood(3,nIPneighbors,Nmembers)) + allocate(geom(ph)%IPareaNormal(3,nIPneighbors,Nmembers)) + allocate(geom(ph)%IParea(nIPneighbors,Nmembers)) call storeGeometry(ph) if(plasticState(ph)%nonlocal .and. .not. allocated(IPneighborhood)) & @@ -652,8 +658,8 @@ module subroutine nonlocal_dependentState(ph, en, ip, el) nRealNeighbors = 0.0_pReal neighbor_rhoTotal = 0.0_pReal do n = 1,nIPneighbors - neighbor_el = IPneighborhood(1,n,ip,el) - neighbor_ip = IPneighborhood(2,n,ip,el) + neighbor_el = geom(ph)%IPneighborhood(1,n,en) + neighbor_ip = geom(ph)%IPneighborhood(2,n,en) no = material_phasememberAt(1,neighbor_ip,neighbor_el) if (neighbor_el > 0 .and. neighbor_ip > 0) then if (material_phaseAt(1,neighbor_el) == ph) then @@ -669,9 +675,9 @@ module subroutine nonlocal_dependentState(ph, en, ip, el) connection_latticeConf(1:3,n) = matmul(invFe, discretization_IPcoords(1:3,neighbor_el+neighbor_ip-1) & - discretization_IPcoords(1:3,el+neighbor_ip-1)) - normal_latticeConf = matmul(transpose(invFp), IPareaNormal(1:3,n,ip,el)) + normal_latticeConf = matmul(transpose(invFp), geom(ph)%IPareaNormal(1:3,n,en)) if (math_inner(normal_latticeConf,connection_latticeConf(1:3,n)) < 0.0_pReal) & ! neighboring connection points in opposite direction to face normal: must be periodic image - connection_latticeConf(1:3,n) = normal_latticeConf * geom(ph)%V_0(en)/IParea(n,ip,el) ! instead take the surface normal scaled with the diameter of the cell + connection_latticeConf(1:3,n) = normal_latticeConf * geom(ph)%V_0(en)/geom(ph)%IParea(n,en) ! instead take the surface normal scaled with the diameter of the cell else ! local neighbor or different lattice structure or different constitution instance -> use central values instead connection_latticeConf(1:3,n) = 0.0_pReal @@ -1101,7 +1107,7 @@ module subroutine nonlocal_dotState(Mp,timestep, & - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) & - rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have - rhoDot = rhoDotFlux(timestep, ph,en,ip,el) & + rhoDot = rhoDotFlux(timestep, ph,en) & + rhoDotMultiplication & + rhoDotSingle2DipoleGlide & + rhoDotAthermalAnnihilation & @@ -1131,17 +1137,15 @@ end subroutine nonlocal_dotState !> @brief calculates the rate of change of microstructure !--------------------------------------------------------------------------------------------------- #if __INTEL_COMPILER >= 2020 -non_recursive function rhoDotFlux(timestep,ph,en,ip,el) +non_recursive function rhoDotFlux(timestep,ph,en) #else -function rhoDotFlux(timestep,ph,en,ip,el) +function rhoDotFlux(timestep,ph,en) #endif real(pReal), intent(in) :: & timestep !< substepped crystallite time increment integer, intent(in) :: & ph, & - en, & - ip, & !< current integration point - el !< current element number + en integer :: & neighbor_ph, & !< phase of my neighbor's plasticity @@ -1217,14 +1221,14 @@ function rhoDotFlux(timestep,ph,en,ip,el) !*** check CFL (Courant-Friedrichs-Lewy) condition for flux if (any( abs(dot_gamma) > 0.0_pReal & ! any active slip system ... .and. prm%C_CFL * abs(v0) * timestep & - > geom(ph)%V_0(en)/ maxval(IParea(:,ip,el)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here) + > geom(ph)%V_0(en)/ maxval(geom(ph)%IParea(:,en)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here) #ifdef DEBUG if (debugConstitutive%extensive) then print'(a,i5,a,i2)', '<< CONST >> CFL condition not fullfilled at ph ',ph,' en ',en print'(a,e10.3,a,e10.3)', '<< CONST >> velocity is at ', & maxval(abs(v0), abs(dot_gamma) > 0.0_pReal & .and. prm%C_CFL * abs(v0) * timestep & - > geom(ph)%V_0(en) / maxval(IParea(:,ip,el))), & + > geom(ph)%V_0(en) / maxval(geom(ph)%IParea(:,en))), & ' at a timestep of ',timestep print*, '<< CONST >> enforcing cutback !!!' end if @@ -1247,16 +1251,16 @@ function rhoDotFlux(timestep,ph,en,ip,el) neighbors: do n = 1,nIPneighbors - neighbor_el = IPneighborhood(1,n,ip,el) - neighbor_ip = IPneighborhood(2,n,ip,el) - neighbor_n = IPneighborhood(3,n,ip,el) + neighbor_el = geom(ph)%IPneighborhood(1,n,en) + neighbor_ip = geom(ph)%IPneighborhood(2,n,en) + neighbor_n = geom(ph)%IPneighborhood(3,n,en) np = material_phaseAt(1,neighbor_el) no = material_phasememberAt(1,neighbor_ip,neighbor_el) opposite_neighbor = n + mod(n,2) - mod(n+1,2) - opposite_el = IPneighborhood(1,opposite_neighbor,ip,el) - opposite_ip = IPneighborhood(2,opposite_neighbor,ip,el) - opposite_n = IPneighborhood(3,opposite_neighbor,ip,el) + opposite_el = geom(ph)%IPneighborhood(1,opposite_neighbor,en) + opposite_ip = geom(ph)%IPneighborhood(2,opposite_neighbor,en) + opposite_n = geom(ph)%IPneighborhood(3,opposite_neighbor,en) if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient neighbor_ph = material_phaseAt(1,neighbor_el) @@ -1327,10 +1331,10 @@ function rhoDotFlux(timestep,ph,en,ip,el) if (phase_plasticity(material_phaseAt(1,opposite_el)) == PLASTIC_NONLOCAL_ID) then normal_me2neighbor_defConf = math_det33(Favg) & - * matmul(math_inv33(transpose(Favg)),IPareaNormal(1:3,n,ip,el)) ! normal of the interface in (average) deformed configuration (pointing en => neighbor) + * matmul(math_inv33(transpose(Favg)),geom(ph)%IPareaNormal(1:3,n,en)) ! normal of the interface in (average) deformed configuration (pointing en => neighbor) normal_me2neighbor = matmul(transpose(my_Fe), normal_me2neighbor_defConf) & / math_det33(my_Fe) ! interface normal in my lattice configuration - area = IParea(n,ip,el) * norm2(normal_me2neighbor) + area = geom(ph)%IParea(n,en) * norm2(normal_me2neighbor) normal_me2neighbor = normal_me2neighbor / norm2(normal_me2neighbor) ! normalize the surface normal to unit length do s = 1,ns do t = 1,4 @@ -1758,14 +1762,27 @@ subroutine storeGeometry(ph) integer, intent(in) :: ph - integer :: ce, co + integer :: ce, co, nCell real(pReal), dimension(:), allocatable :: V + integer, dimension(:,:,:), allocatable :: neighborhood + real(pReal), dimension(:,:), allocatable :: area + real(pReal), dimension(:,:,:), allocatable :: areaNormal + nCell = product(shape(IPvolume)) + + V = reshape(IPvolume,[nCell]) + neighborhood = reshape(IPneighborhood,[3,nIPneighbors,nCell]) + area = reshape(IParea,[nIPneighbors,nCell]) + areaNormal = reshape(IPareaNormal,[3,nIPneighbors,nCell]) - V = reshape(IPvolume,[product(shape(IPvolume))]) do ce = 1, size(material_homogenizationEntry,1) do co = 1, homogenization_maxNconstituents - if (material_phaseID(co,ce) == ph) geom(ph)%V_0(material_phaseEntry(co,ce)) = V(ce) + if (material_phaseID(co,ce) == ph) then + geom(ph)%V_0(material_phaseEntry(co,ce)) = V(ce) + geom(ph)%IPneighborhood(:,:,material_phaseEntry(co,ce)) = neighborhood(:,:,ce) + geom(ph)%IParea(:,material_phaseEntry(co,ce)) = area(:,ce) + geom(ph)%IPareaNormal(:,:,material_phaseEntry(co,ce)) = areaNormal(:,:,ce) + end if end do end do From 2bd10a126178a540b7c85e8f72e08fe80c674537 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 4 Feb 2022 08:25:51 +0100 Subject: [PATCH 04/13] avoid use of IP/el --- src/phase_mechanical_plastic.f90 | 16 ++++++---------- src/phase_mechanical_plastic_nonlocal.f90 | 23 +++++++++++------------ 2 files changed, 17 insertions(+), 22 deletions(-) diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index 912aadc03..9378edfdc 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -154,16 +154,14 @@ submodule(phase:mechanical) plastic en end subroutine dislotungsten_dotState - module subroutine nonlocal_dotState(Mp,timestep,ph,en,ip,el) + module subroutine nonlocal_dotState(Mp,timestep,ph,en) real(pReal), dimension(3,3), intent(in) :: & Mp !< MandelStress real(pReal), intent(in) :: & timestep !< substepped crystallite time increment integer, intent(in) :: & ph, & - en, & - ip, & !< current integration point - el !< current element number + en end subroutine nonlocal_dotState module subroutine dislotwin_dependentState(T,ph,en) @@ -180,12 +178,10 @@ submodule(phase:mechanical) plastic en end subroutine dislotungsten_dependentState - module subroutine nonlocal_dependentState(ph, en, ip, el) + module subroutine nonlocal_dependentState(ph, en) integer, intent(in) :: & ph, & - en, & - ip, & !< current integration point - el !< current element number + en end subroutine nonlocal_dependentState module subroutine plastic_kinehardening_deltaState(Mp,ph,en) @@ -333,7 +329,7 @@ module function plastic_dotState(subdt,co,ip,el,ph,en) result(dotState) call dislotungsten_dotState(Mp,thermal_T(ph,en),ph,en) case (PLASTIC_NONLOCAL_ID) plasticType - call nonlocal_dotState(Mp,subdt,ph,en,ip,el) + call nonlocal_dotState(Mp,subdt,ph,en) end select plasticType end if @@ -369,7 +365,7 @@ module subroutine plastic_dependentState(co, ip, el) call dislotungsten_dependentState(ph,en) case (PLASTIC_NONLOCAL_ID) plasticType - call nonlocal_dependentState(ph,en,ip,el) + call nonlocal_dependentState(ph,en) end select plasticType diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 78271b92f..9c901f6b0 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -16,7 +16,7 @@ submodule(phase:plastic) nonlocal type :: tGeometry real(pReal), dimension(:), allocatable :: V_0 integer, dimension(:,:,:), allocatable :: IPneighborhood - real(pReal), dimension(:,:), allocatable :: IParea + real(pReal), dimension(:,:), allocatable :: IParea, IPcoordinates real(pReal), dimension(:,:,:), allocatable :: IPareaNormal end type tGeometry @@ -413,6 +413,7 @@ module function plastic_nonlocal_init() result(myPlasticity) allocate(geom(ph)%IPneighborhood(3,nIPneighbors,Nmembers)) allocate(geom(ph)%IPareaNormal(3,nIPneighbors,Nmembers)) allocate(geom(ph)%IParea(nIPneighbors,Nmembers)) + allocate(geom(ph)%IPcoordinates(3,Nmembers)) call storeGeometry(ph) if(plasticState(ph)%nonlocal .and. .not. allocated(IPneighborhood)) & @@ -551,13 +552,11 @@ end function plastic_nonlocal_init !-------------------------------------------------------------------------------------------------- !> @brief calculates quantities characterizing the microstructure !-------------------------------------------------------------------------------------------------- -module subroutine nonlocal_dependentState(ph, en, ip, el) +module subroutine nonlocal_dependentState(ph, en) integer, intent(in) :: & ph, & - en, & - ip, & - el + en integer :: & no, & !< neighbor offset @@ -673,8 +672,8 @@ module subroutine nonlocal_dependentState(ph, en, ip, el) neighbor_rhoTotal(1,:,n) = sum(abs(rho_neighbor0(:,edg)),2) neighbor_rhoTotal(2,:,n) = sum(abs(rho_neighbor0(:,scr)),2) - connection_latticeConf(1:3,n) = matmul(invFe, discretization_IPcoords(1:3,neighbor_el+neighbor_ip-1) & - - discretization_IPcoords(1:3,el+neighbor_ip-1)) + connection_latticeConf(1:3,n) = matmul(invFe, geom(ph)%IPcoordinates(1:3,en) & + - geom(ph)%IPcoordinates(1:3,en-1)) ! ToDo: broken for different materials normal_latticeConf = matmul(transpose(invFp), geom(ph)%IPareaNormal(1:3,n,en)) if (math_inner(normal_latticeConf,connection_latticeConf(1:3,n)) < 0.0_pReal) & ! neighboring connection points in opposite direction to face normal: must be periodic image connection_latticeConf(1:3,n) = normal_latticeConf * geom(ph)%V_0(en)/geom(ph)%IParea(n,en) ! instead take the surface normal scaled with the diameter of the cell @@ -947,7 +946,7 @@ end subroutine plastic_nonlocal_deltaState !> @brief calculates the rate of change of microstructure !--------------------------------------------------------------------------------------------------- module subroutine nonlocal_dotState(Mp,timestep, & - ph,en,ip,el) + ph,en) real(pReal), dimension(3,3), intent(in) :: & Mp !< MandelStress @@ -955,9 +954,7 @@ module subroutine nonlocal_dotState(Mp,timestep, & timestep !< substepped crystallite time increment integer, intent(in) :: & ph, & - en, & - ip, & !< current integration point - el !< current element number + en integer :: & c, & !< character of dislocation @@ -1765,7 +1762,7 @@ subroutine storeGeometry(ph) integer :: ce, co, nCell real(pReal), dimension(:), allocatable :: V integer, dimension(:,:,:), allocatable :: neighborhood - real(pReal), dimension(:,:), allocatable :: area + real(pReal), dimension(:,:), allocatable :: area, coords real(pReal), dimension(:,:,:), allocatable :: areaNormal nCell = product(shape(IPvolume)) @@ -1774,6 +1771,7 @@ subroutine storeGeometry(ph) neighborhood = reshape(IPneighborhood,[3,nIPneighbors,nCell]) area = reshape(IParea,[nIPneighbors,nCell]) areaNormal = reshape(IPareaNormal,[3,nIPneighbors,nCell]) + coords = reshape(discretization_IPcoords,[3,nCell]) do ce = 1, size(material_homogenizationEntry,1) do co = 1, homogenization_maxNconstituents @@ -1782,6 +1780,7 @@ subroutine storeGeometry(ph) geom(ph)%IPneighborhood(:,:,material_phaseEntry(co,ce)) = neighborhood(:,:,ce) geom(ph)%IParea(:,material_phaseEntry(co,ce)) = area(:,ce) geom(ph)%IPareaNormal(:,:,material_phaseEntry(co,ce)) = areaNormal(:,:,ce) + geom(ph)%IPcoordinates(:,material_phaseEntry(co,ce)) = coords(:,ce) end if end do end do From 3a0596a274f673a3bf93b308c6b4cf2b41c4835f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 4 Feb 2022 08:49:00 +0100 Subject: [PATCH 05/13] use en/ph access --- src/phase.f90 | 15 ++++++++------- src/phase_mechanical.f90 | 2 +- src/phase_mechanical_plastic.f90 | 14 +++----------- src/phase_mechanical_plastic_nonlocal.f90 | 4 ++-- 4 files changed, 14 insertions(+), 21 deletions(-) diff --git a/src/phase.f90 b/src/phase.f90 index 82bc44bfd..4d180d439 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -272,14 +272,12 @@ module phase type(tRotationContainer), dimension(:), intent(in) :: orientation end subroutine plastic_nonlocal_updateCompatibility - module subroutine plastic_dependentState(co,ip,el) + module subroutine plastic_dependentState(en,ph) integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element + en, & + ph end subroutine plastic_dependentState - module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,en) integer, intent(in) :: ph, en real(pReal), intent(in), dimension(3,3) :: & @@ -503,7 +501,8 @@ subroutine crystallite_init() el, & !< counter in element loop cMax, & !< maximum number of integration point components iMax, & !< maximum number of integration points - eMax !< maximum number of elements + eMax, & !< maximum number of elements + en, ph class(tNode), pointer :: & num_crystallite, & @@ -560,8 +559,10 @@ subroutine crystallite_init() do ip = 1, iMax ce = (el-1)*discretization_nIPs + ip do co = 1,homogenization_Nconstituents(material_homogenizationID(ce)) + en = material_phaseEntry(co,ce) + ph = material_phaseID(co,ce) call crystallite_orientations(co,ip,el) - call plastic_dependentState(co,ip,el) ! update dependent state variables to be consistent with basic states + call plastic_dependentState(en,ph) ! update dependent state variables to be consistent with basic states end do end do end do diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 04cc4b946..c97e4a13b 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -431,7 +431,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken) ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) - call plastic_dependentState(co,ip,el) + call plastic_dependentState(en,ph) Lpguess = phase_mechanical_Lp(ph)%data(1:3,1:3,en) ! take as first guess Liguess = phase_mechanical_Li(ph)%data(1:3,1:3,en) ! take as first guess diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index 9378edfdc..a5f40a2df 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -341,20 +341,12 @@ end function plastic_dotState !-------------------------------------------------------------------------------------------------- !> @brief calls microstructure function of the different plasticity constitutive models !-------------------------------------------------------------------------------------------------- -module subroutine plastic_dependentState(co, ip, el) +module subroutine plastic_dependentState(en,ph) integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element + en, & + ph - integer :: & - ph, & - en - - - ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) - en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) plasticType: select case (phase_plasticity(ph)) diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 9c901f6b0..dec7ea8e7 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -672,8 +672,8 @@ module subroutine nonlocal_dependentState(ph, en) neighbor_rhoTotal(1,:,n) = sum(abs(rho_neighbor0(:,edg)),2) neighbor_rhoTotal(2,:,n) = sum(abs(rho_neighbor0(:,scr)),2) - connection_latticeConf(1:3,n) = matmul(invFe, geom(ph)%IPcoordinates(1:3,en) & - - geom(ph)%IPcoordinates(1:3,en-1)) ! ToDo: broken for different materials + connection_latticeConf(1:3,n) = matmul(invFe, geom(ph)%IPcoordinates(1:3,no) & + - geom(ph)%IPcoordinates(1:3,en)) normal_latticeConf = matmul(transpose(invFp), geom(ph)%IPareaNormal(1:3,n,en)) if (math_inner(normal_latticeConf,connection_latticeConf(1:3,n)) < 0.0_pReal) & ! neighboring connection points in opposite direction to face normal: must be periodic image connection_latticeConf(1:3,n) = normal_latticeConf * geom(ph)%V_0(en)/geom(ph)%IParea(n,en) ! instead take the surface normal scaled with the diameter of the cell From 8aed927989ee5edec91f03479bb110e660c3d4c6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 4 Feb 2022 11:28:54 +0100 Subject: [PATCH 06/13] avoid use of co/ip/el at the phase level --- src/phase_mechanical.f90 | 105 ++++++++++++------------------- src/phase_mechanical_plastic.f90 | 5 +- 2 files changed, 42 insertions(+), 68 deletions(-) diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index c97e4a13b..88e98cc2e 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -79,15 +79,12 @@ submodule(phase) mechanical en end subroutine plastic_isotropic_LiAndItsTangent - module function plastic_dotState(subdt,co,ip,el,ph,en) result(dotState) + module function plastic_dotState(subdt,ph,en) result(dotState) integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el, & !< element ph, & en real(pReal), intent(in) :: & - subdt !< timestep + subdt !< timestep real(pReal), dimension(plasticState(ph)%sizeDotState) :: & dotState end function plastic_dotState @@ -365,13 +362,11 @@ end subroutine mechanical_results !> @brief calculation of stress (P) with time integration based on a residuum in Lp and !> intermediate acceleration of the Newton-Raphson correction !-------------------------------------------------------------------------------------------------- -function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken) +function integrateStress(F,subFp0,subFi0,Delta_t,en,ph) result(broken) real(pReal), dimension(3,3), intent(in) :: F,subFp0,subFi0 real(pReal), intent(in) :: Delta_t - integer, intent(in):: el, & ! element index - ip, & ! integration point index - co ! grain index + integer, intent(in) :: en, ph real(pReal), dimension(3,3):: Fp_new, & ! plastic deformation gradient at end of timestep invFp_new, & ! inverse of Fp_new @@ -427,10 +422,6 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken) broken = .true. - - ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) - en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) - call plastic_dependentState(en,ph) Lpguess = phase_mechanical_Lp(ph)%data(1:3,1:3,en) ! take as first guess @@ -579,15 +570,14 @@ end function integrateStress !> @brief integrate stress, state with adaptive 1st order explicit Euler method !> using Fixed Point Iteration to adapt the stepsize !-------------------------------------------------------------------------------------------------- -function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) result(broken) +function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,en,ph) result(broken) real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0 real(pReal), intent(in),dimension(:) :: subState0 real(pReal), intent(in) :: Delta_t integer, intent(in) :: & - el, & !< element index in element loop - ip, & !< integration point index in ip loop - co !< grain index in grain loop + en, & + ph logical :: & broken @@ -598,18 +588,16 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul sizeDotState real(pReal) :: & zeta - real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: & + real(pReal), dimension(plasticState(ph)%sizeDotState) :: & r, & ! state residuum dotState - real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState,2) :: & + real(pReal), dimension(plasticState(ph)%sizeDotState,2) :: & dotState_last - ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) - en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) broken = .true. - dotState = plastic_dotState(Delta_t, co,ip,el,ph,en) + dotState = plastic_dotState(Delta_t,ph,en) if (any(IEEE_is_NaN(dotState))) return sizeDotState = plasticState(ph)%sizeDotState @@ -620,10 +608,10 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul dotState_last(1:sizeDotState,2) = merge(dotState_last(1:sizeDotState,1),0.0, nIterationState > 1) dotState_last(1:sizeDotState,1) = dotState - broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) + broken = integrateStress(F,subFp0,subFi0,Delta_t,en,ph) if(broken) exit iteration - dotState = plastic_dotState(Delta_t, co,ip,el,ph,en) + dotState = plastic_dotState(Delta_t,ph,en) if (any(IEEE_is_NaN(dotState))) exit iteration zeta = damper(dotState,dotState_last(1:sizeDotState,1),dotState_last(1:sizeDotState,2)) @@ -672,19 +660,18 @@ end function integrateStateFPI !-------------------------------------------------------------------------------------------------- !> @brief integrate state with 1st order explicit Euler method !-------------------------------------------------------------------------------------------------- -function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) result(broken) +function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,en,ph) result(broken) real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0 real(pReal), intent(in),dimension(:) :: subState0 real(pReal), intent(in) :: Delta_t integer, intent(in) :: & - el, & !< element index in element loop - ip, & !< integration point index in ip loop - co !< grain index in grain loop + en, & + ph !< grain index in grain loop logical :: & broken - real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: & + real(pReal), dimension(plasticState(ph)%sizeDotState) :: & dotState integer :: & ph, & @@ -692,11 +679,9 @@ function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) res sizeDotState - ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) - en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) broken = .true. - dotState = plastic_dotState(Delta_t, co,ip,el,ph,en) + dotState = plastic_dotState(Delta_t,ph,en) if (any(IEEE_is_NaN(dotState))) return sizeDotState = plasticState(ph)%sizeDotState @@ -706,7 +691,7 @@ function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) res broken = plastic_deltaState(ph,en) if(broken) return - broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) + broken = integrateStress(F,subFp0,subFi0,Delta_t,en,ph) end function integrateStateEuler @@ -714,15 +699,14 @@ end function integrateStateEuler !-------------------------------------------------------------------------------------------------- !> @brief integrate stress, state with 1st order Euler method with adaptive step size !-------------------------------------------------------------------------------------------------- -function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) result(broken) +function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,en,ph) result(broken) real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0 real(pReal), intent(in),dimension(:) :: subState0 real(pReal), intent(in) :: Delta_t integer, intent(in) :: & - el, & !< element index in element loop - ip, & !< integration point index in ip loop - co !< grain index in grain loop + en, & + ph logical :: & broken @@ -730,16 +714,14 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip ph, & en, & sizeDotState - real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: & + real(pReal), dimension(plasticState(ph)%sizeDotState) :: & r, & dotState - ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) - en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) broken = .true. - dotState = plastic_dotState(Delta_t, co,ip,el,ph,en) + dotState = plastic_dotState(Delta_t,ph,en) if (any(IEEE_is_NaN(dotState))) return sizeDotState = plasticState(ph)%sizeDotState @@ -751,10 +733,10 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip broken = plastic_deltaState(ph,en) if(broken) return - broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) + broken = integrateStress(F,subFp0,subFi0,Delta_t,en,ph) if(broken) return - dotState = plastic_dotState(Delta_t, co,ip,el,ph,en) + dotState = plastic_dotState(Delta_t,ph,en) if (any(IEEE_is_NaN(dotState))) return broken = .not. converged(r + 0.5_pReal * dotState * Delta_t, & @@ -767,12 +749,12 @@ end function integrateStateAdaptiveEuler !--------------------------------------------------------------------------------------------------- !> @brief Integrate state (including stress integration) with the classic Runge Kutta method !--------------------------------------------------------------------------------------------------- -function integrateStateRK4(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) result(broken) +function integrateStateRK4(F_0,F,subFp0,subFi0,subState0,Delta_t,en,ph) result(broken) real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0 real(pReal), intent(in),dimension(:) :: subState0 real(pReal), intent(in) :: Delta_t - integer, intent(in) :: co,ip,el + integer, intent(in) :: en, ph logical :: broken real(pReal), dimension(3,3), parameter :: & @@ -787,7 +769,7 @@ function integrateStateRK4(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul B = [1.0_pReal/6.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/6.0_pReal] - broken = integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C) + broken = integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,en,ph,A,B,C) end function integrateStateRK4 @@ -795,12 +777,12 @@ end function integrateStateRK4 !--------------------------------------------------------------------------------------------------- !> @brief Integrate state (including stress integration) with the Cash-Carp method !--------------------------------------------------------------------------------------------------- -function integrateStateRKCK45(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) result(broken) +function integrateStateRKCK45(F_0,F,subFp0,subFi0,subState0,Delta_t,en,ph) result(broken) real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0 real(pReal), intent(in),dimension(:) :: subState0 real(pReal), intent(in) :: Delta_t - integer, intent(in) :: co,ip,el + integer, intent(in) :: en, ph logical :: broken real(pReal), dimension(5,5), parameter :: & @@ -822,7 +804,7 @@ function integrateStateRKCK45(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) re 13525.0_pReal/55296.0_pReal, 277.0_pReal/14336.0_pReal, 1._pReal/4._pReal] - broken = integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,DB) + broken = integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,en,ph,A,B,C,DB) end function integrateStateRKCK45 @@ -831,7 +813,7 @@ end function integrateStateRKCK45 !> @brief Integrate state (including stress integration) with an explicit Runge-Kutta method or an !! embedded explicit Runge-Kutta method !-------------------------------------------------------------------------------------------------- -function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,DB) result(broken) +function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,en,ph,A,B,C,DB) result(broken) real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0 real(pReal), intent(in),dimension(:) :: subState0 @@ -840,9 +822,8 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D real(pReal), dimension(:), intent(in) :: B, C real(pReal), dimension(:), intent(in), optional :: DB integer, intent(in) :: & - el, & !< element index in element loop - ip, & !< integration point index in ip loop - co !< grain index in grain loop + en, & + ph logical :: broken integer :: & @@ -851,17 +832,15 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D ph, & en, & sizeDotState - real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: & + real(pReal), dimension(plasticState(ph)%sizeDotState) :: & dotState - real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState,size(B)) :: & + real(pReal), dimension(plasticState(ph)%sizeDotState,size(B)) :: & plastic_RKdotState - ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) - en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) broken = .true. - dotState = plastic_dotState(Delta_t, co,ip,el,ph,en) + dotState = plastic_dotState(Delta_t,ph,en) if (any(IEEE_is_NaN(dotState))) return sizeDotState = plasticState(ph)%sizeDotState @@ -879,10 +858,10 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D plasticState(ph)%state(1:sizeDotState,en) = subState0 & + dotState * Delta_t - broken = integrateStress(F_0 + (F - F_0) * Delta_t * C(stage),subFp0,subFi0,Delta_t * C(stage),co,ip,el) + broken = integrateStress(F_0 + (F - F_0) * Delta_t * C(stage),subFp0,subFi0,Delta_t * C(stage),en,ph) if(broken) exit - dotState = plastic_dotState(Delta_t, co,ip,el,ph,en) + dotState = plastic_dotState(Delta_t,ph,en) if (any(IEEE_is_NaN(dotState))) exit enddo @@ -904,7 +883,7 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D broken = plastic_deltaState(ph,en) if(broken) return - broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) + broken = integrateStress(F,subFp0,subFi0,Delta_t,en,ph) end function integrateStateRK @@ -1031,8 +1010,6 @@ module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged real(pReal), dimension(:), allocatable :: subState0 - ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) - en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) sizeDotState = plasticState(ph)%sizeDotState subLi0 = phase_mechanical_Li0(ph)%data(1:3,1:3,en) @@ -1084,7 +1061,7 @@ module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged if (todo) then subF = subF0 & + subStep * (phase_mechanical_F(ph)%data(1:3,1:3,en) - phase_mechanical_F0(ph)%data(1:3,1:3,en)) - converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * Delta_t,co,ip,el) + converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * Delta_t,en,ph) endif enddo cutbackLooping diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index a5f40a2df..13f14717e 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -291,12 +291,9 @@ end subroutine plastic_LpAndItsTangents !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -module function plastic_dotState(subdt,co,ip,el,ph,en) result(dotState) +module function plastic_dotState(subdt,ph,en) result(dotState) integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el, & !< element ph, & en real(pReal), intent(in) :: & From afed13e2ca5d6c33bc5b85e339512b99e29eca03 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 4 Feb 2022 12:25:11 +0100 Subject: [PATCH 07/13] ip/el not of interest --- src/homogenization.f90 | 4 ++-- src/phase.f90 | 8 ++++---- src/phase_damage.f90 | 9 ++++----- src/phase_mechanical.f90 | 17 ++++------------- 4 files changed, 14 insertions(+), 24 deletions(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 3d96b007f..1f9fb2221 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -259,7 +259,7 @@ subroutine homogenization_mechanical_response(Delta_t,FEsolving_execIP,FEsolving NiterationMPstate = NiterationMPstate + 1 call mechanical_partition(homogenization_F(1:3,1:3,ce),ce) - converged = all([(phase_mechanical_constitutive(Delta_t,co,ip,el),co=1,homogenization_Nconstituents(ho))]) + converged = all([(phase_mechanical_constitutive(Delta_t,co,ce),co=1,homogenization_Nconstituents(ho))]) if (converged) then doneAndHappy = mechanical_updateState(Delta_t,homogenization_F(1:3,1:3,ce),ce) converged = all(doneAndHappy) @@ -268,7 +268,7 @@ subroutine homogenization_mechanical_response(Delta_t,FEsolving_execIP,FEsolving endif enddo convergenceLooping - converged = converged .and. all([(phase_damage_constitutive(Delta_t,co,ip,el),co=1,homogenization_Nconstituents(ho))]) + converged = converged .and. all([(phase_damage_constitutive(Delta_t,co,ce),co=1,homogenization_Nconstituents(ho))]) if (.not. converged) then if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill' diff --git a/src/phase.f90 b/src/phase.f90 index 4d180d439..8bc797e85 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -228,15 +228,15 @@ module phase end function phase_thermal_constitutive - module function phase_damage_constitutive(Delta_t,co,ip,el) result(converged_) + module function phase_damage_constitutive(Delta_t,co,ce) result(converged_) real(pReal), intent(in) :: Delta_t - integer, intent(in) :: co, ip, el + integer, intent(in) :: co, ce logical :: converged_ end function phase_damage_constitutive - module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged_) + module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_) real(pReal), intent(in) :: Delta_t - integer, intent(in) :: co, ip, el + integer, intent(in) :: co, ce logical :: converged_ end function phase_mechanical_constitutive diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 2fa5ab9ea..e4700938f 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -127,21 +127,20 @@ end subroutine damage_init !-------------------------------------------------------------------------------------------------- !> @brief calculate stress (P) !-------------------------------------------------------------------------------------------------- -module function phase_damage_constitutive(Delta_t,co,ip,el) result(converged_) +module function phase_damage_constitutive(Delta_t,co,ce) result(converged_) real(pReal), intent(in) :: Delta_t integer, intent(in) :: & co, & - ip, & - el + ce logical :: converged_ integer :: & ph, en - ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) - en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) + ph = material_phaseID(co,ce) + en = material_phaseEntry(co,ce) converged_ = .not. integrateDamageState(Delta_t,ph,en) diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 88e98cc2e..655f53ad8 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -414,8 +414,6 @@ function integrateStress(F,subFp0,subFi0,Delta_t,en,ph) result(broken) ierr, & ! error indicator for LAPACK o, & p, & - ph, & - en, & jacoCounterLp, & jacoCounterLi ! counters to check for Jacobian update logical :: error,broken @@ -583,8 +581,6 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,en,ph) result(b integer :: & NiterationState, & !< number of iterations in state loop - ph, & - en, & sizeDotState real(pReal) :: & zeta @@ -674,8 +670,6 @@ function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,en,ph) result real(pReal), dimension(plasticState(ph)%sizeDotState) :: & dotState integer :: & - ph, & - en, & sizeDotState @@ -711,8 +705,6 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,en,ph broken integer :: & - ph, & - en, & sizeDotState real(pReal), dimension(plasticState(ph)%sizeDotState) :: & r, & @@ -829,8 +821,6 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,en,ph,A,B,C,DB) integer :: & stage, & ! stage index in integration stage loop n, & - ph, & - en, & sizeDotState real(pReal), dimension(plasticState(ph)%sizeDotState) :: & dotState @@ -985,13 +975,12 @@ end subroutine mechanical_forward !-------------------------------------------------------------------------------------------------- !> @brief calculate stress (P) !-------------------------------------------------------------------------------------------------- -module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged_) +module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_) real(pReal), intent(in) :: Delta_t integer, intent(in) :: & co, & - ip, & - el + ce logical :: converged_ real(pReal) :: & @@ -1010,6 +999,8 @@ module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged real(pReal), dimension(:), allocatable :: subState0 + ph = material_phaseID(co,ce) + en = material_phaseEntry(co,ce) sizeDotState = plasticState(ph)%sizeDotState subLi0 = phase_mechanical_Li0(ph)%data(1:3,1:3,en) From 3990536f1ca90d6e2f9709bd8b9b2f54732819e8 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 4 Feb 2022 17:42:05 +0100 Subject: [PATCH 08/13] consistent order of arguments --- src/phase.f90 | 8 ++--- src/phase_mechanical.f90 | 54 ++++++++++++++++---------------- src/phase_mechanical_plastic.f90 | 8 ++--- 3 files changed, 35 insertions(+), 35 deletions(-) diff --git a/src/phase.f90 b/src/phase.f90 index 8bc797e85..c4e7142f8 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -272,10 +272,10 @@ module phase type(tRotationContainer), dimension(:), intent(in) :: orientation end subroutine plastic_nonlocal_updateCompatibility - module subroutine plastic_dependentState(en,ph) + module subroutine plastic_dependentState(ph,en) integer, intent(in) :: & - en, & - ph + ph, & + en end subroutine plastic_dependentState module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,en) @@ -562,7 +562,7 @@ subroutine crystallite_init() en = material_phaseEntry(co,ce) ph = material_phaseID(co,ce) call crystallite_orientations(co,ip,el) - call plastic_dependentState(en,ph) ! update dependent state variables to be consistent with basic states + call plastic_dependentState(ph,en) ! update dependent state variables to be consistent with basic states end do end do end do diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 655f53ad8..de9c2bea2 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -362,11 +362,11 @@ end subroutine mechanical_results !> @brief calculation of stress (P) with time integration based on a residuum in Lp and !> intermediate acceleration of the Newton-Raphson correction !-------------------------------------------------------------------------------------------------- -function integrateStress(F,subFp0,subFi0,Delta_t,en,ph) result(broken) +function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken) real(pReal), dimension(3,3), intent(in) :: F,subFp0,subFi0 real(pReal), intent(in) :: Delta_t - integer, intent(in) :: en, ph + integer, intent(in) :: ph, en real(pReal), dimension(3,3):: Fp_new, & ! plastic deformation gradient at end of timestep invFp_new, & ! inverse of Fp_new @@ -420,7 +420,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,en,ph) result(broken) broken = .true. - call plastic_dependentState(en,ph) + call plastic_dependentState(ph,en) Lpguess = phase_mechanical_Lp(ph)%data(1:3,1:3,en) ! take as first guess Liguess = phase_mechanical_Li(ph)%data(1:3,1:3,en) ! take as first guess @@ -568,14 +568,14 @@ end function integrateStress !> @brief integrate stress, state with adaptive 1st order explicit Euler method !> using Fixed Point Iteration to adapt the stepsize !-------------------------------------------------------------------------------------------------- -function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,en,ph) result(broken) +function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(broken) real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0 real(pReal), intent(in),dimension(:) :: subState0 real(pReal), intent(in) :: Delta_t integer, intent(in) :: & - en, & - ph + ph, & + en logical :: & broken @@ -604,7 +604,7 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,en,ph) result(b dotState_last(1:sizeDotState,2) = merge(dotState_last(1:sizeDotState,1),0.0, nIterationState > 1) dotState_last(1:sizeDotState,1) = dotState - broken = integrateStress(F,subFp0,subFi0,Delta_t,en,ph) + broken = integrateStress(F,subFp0,subFi0,Delta_t,ph,en) if(broken) exit iteration dotState = plastic_dotState(Delta_t,ph,en) @@ -656,14 +656,14 @@ end function integrateStateFPI !-------------------------------------------------------------------------------------------------- !> @brief integrate state with 1st order explicit Euler method !-------------------------------------------------------------------------------------------------- -function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,en,ph) result(broken) +function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(broken) real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0 real(pReal), intent(in),dimension(:) :: subState0 real(pReal), intent(in) :: Delta_t integer, intent(in) :: & - en, & - ph !< grain index in grain loop + ph, & + en !< grain index in grain loop logical :: & broken @@ -685,7 +685,7 @@ function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,en,ph) result broken = plastic_deltaState(ph,en) if(broken) return - broken = integrateStress(F,subFp0,subFi0,Delta_t,en,ph) + broken = integrateStress(F,subFp0,subFi0,Delta_t,ph,en) end function integrateStateEuler @@ -693,14 +693,14 @@ end function integrateStateEuler !-------------------------------------------------------------------------------------------------- !> @brief integrate stress, state with 1st order Euler method with adaptive step size !-------------------------------------------------------------------------------------------------- -function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,en,ph) result(broken) +function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(broken) real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0 real(pReal), intent(in),dimension(:) :: subState0 real(pReal), intent(in) :: Delta_t integer, intent(in) :: & - en, & - ph + ph, & + en logical :: & broken @@ -725,7 +725,7 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,en,ph broken = plastic_deltaState(ph,en) if(broken) return - broken = integrateStress(F,subFp0,subFi0,Delta_t,en,ph) + broken = integrateStress(F,subFp0,subFi0,Delta_t,ph,en) if(broken) return dotState = plastic_dotState(Delta_t,ph,en) @@ -741,12 +741,12 @@ end function integrateStateAdaptiveEuler !--------------------------------------------------------------------------------------------------- !> @brief Integrate state (including stress integration) with the classic Runge Kutta method !--------------------------------------------------------------------------------------------------- -function integrateStateRK4(F_0,F,subFp0,subFi0,subState0,Delta_t,en,ph) result(broken) +function integrateStateRK4(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(broken) real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0 real(pReal), intent(in),dimension(:) :: subState0 real(pReal), intent(in) :: Delta_t - integer, intent(in) :: en, ph + integer, intent(in) :: ph, en logical :: broken real(pReal), dimension(3,3), parameter :: & @@ -761,7 +761,7 @@ function integrateStateRK4(F_0,F,subFp0,subFi0,subState0,Delta_t,en,ph) result(b B = [1.0_pReal/6.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/6.0_pReal] - broken = integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,en,ph,A,B,C) + broken = integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en,A,B,C) end function integrateStateRK4 @@ -769,12 +769,12 @@ end function integrateStateRK4 !--------------------------------------------------------------------------------------------------- !> @brief Integrate state (including stress integration) with the Cash-Carp method !--------------------------------------------------------------------------------------------------- -function integrateStateRKCK45(F_0,F,subFp0,subFi0,subState0,Delta_t,en,ph) result(broken) +function integrateStateRKCK45(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(broken) real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0 real(pReal), intent(in),dimension(:) :: subState0 real(pReal), intent(in) :: Delta_t - integer, intent(in) :: en, ph + integer, intent(in) :: ph, en logical :: broken real(pReal), dimension(5,5), parameter :: & @@ -796,7 +796,7 @@ function integrateStateRKCK45(F_0,F,subFp0,subFi0,subState0,Delta_t,en,ph) resul 13525.0_pReal/55296.0_pReal, 277.0_pReal/14336.0_pReal, 1._pReal/4._pReal] - broken = integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,en,ph,A,B,C,DB) + broken = integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en,A,B,C,DB) end function integrateStateRKCK45 @@ -805,7 +805,7 @@ end function integrateStateRKCK45 !> @brief Integrate state (including stress integration) with an explicit Runge-Kutta method or an !! embedded explicit Runge-Kutta method !-------------------------------------------------------------------------------------------------- -function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,en,ph,A,B,C,DB) result(broken) +function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en,A,B,C,DB) result(broken) real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0 real(pReal), intent(in),dimension(:) :: subState0 @@ -814,8 +814,8 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,en,ph,A,B,C,DB) real(pReal), dimension(:), intent(in) :: B, C real(pReal), dimension(:), intent(in), optional :: DB integer, intent(in) :: & - en, & - ph + ph, & + en logical :: broken integer :: & @@ -848,7 +848,7 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,en,ph,A,B,C,DB) plasticState(ph)%state(1:sizeDotState,en) = subState0 & + dotState * Delta_t - broken = integrateStress(F_0 + (F - F_0) * Delta_t * C(stage),subFp0,subFi0,Delta_t * C(stage),en,ph) + broken = integrateStress(F_0 + (F - F_0) * Delta_t * C(stage),subFp0,subFi0,Delta_t * C(stage),ph,en) if(broken) exit dotState = plastic_dotState(Delta_t,ph,en) @@ -873,7 +873,7 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,en,ph,A,B,C,DB) broken = plastic_deltaState(ph,en) if(broken) return - broken = integrateStress(F,subFp0,subFi0,Delta_t,en,ph) + broken = integrateStress(F,subFp0,subFi0,Delta_t,ph,en) end function integrateStateRK @@ -1052,7 +1052,7 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_) if (todo) then subF = subF0 & + subStep * (phase_mechanical_F(ph)%data(1:3,1:3,en) - phase_mechanical_F0(ph)%data(1:3,1:3,en)) - converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * Delta_t,en,ph) + converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * Delta_t,ph,en) endif enddo cutbackLooping diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index 13f14717e..973f0ddbf 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -178,7 +178,7 @@ submodule(phase:mechanical) plastic en end subroutine dislotungsten_dependentState - module subroutine nonlocal_dependentState(ph, en) + module subroutine nonlocal_dependentState(ph,en) integer, intent(in) :: & ph, & en @@ -338,11 +338,11 @@ end function plastic_dotState !-------------------------------------------------------------------------------------------------- !> @brief calls microstructure function of the different plasticity constitutive models !-------------------------------------------------------------------------------------------------- -module subroutine plastic_dependentState(en,ph) +module subroutine plastic_dependentState(ph,en) integer, intent(in) :: & - en, & - ph + ph, & + en plasticType: select case (phase_plasticity(ph)) From 3c148b5b0e0e4229cafe97adea42d5731ab9056d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 4 Feb 2022 18:55:59 +0100 Subject: [PATCH 09/13] bugfix: openMP variable was not protected --- src/homogenization.f90 | 2 +- src/phase.f90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 1f9fb2221..8c9ecaecb 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -237,7 +237,7 @@ subroutine homogenization_mechanical_response(Delta_t,FEsolving_execIP,FEsolving doneAndHappy - !$OMP PARALLEL DO PRIVATE(ce,en,ho,NiterationMPstate,converged,doneAndHappy) + !$OMP PARALLEL DO PRIVATE(ce,co,en,ho,NiterationMPstate,converged,doneAndHappy) do el = FEsolving_execElem(1),FEsolving_execElem(2) do ip = FEsolving_execIP(1),FEsolving_execIP(2) diff --git a/src/phase.f90 b/src/phase.f90 index c4e7142f8..9943c777a 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -554,7 +554,7 @@ subroutine crystallite_init() flush(IO_STDOUT) - !$OMP PARALLEL DO PRIVATE(ce) + !$OMP PARALLEL DO PRIVATE(ce,ph,en) do el = 1, eMax do ip = 1, iMax ce = (el-1)*discretization_nIPs + ip From 6d78400f875b7642baf9d3ad5ad86111c403ea5e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 5 Feb 2022 07:29:00 +0100 Subject: [PATCH 10/13] the concept of IP/element_ID should not be used at the DAMASK core --- src/CPFEM.f90 | 3 +- src/grid/spectral_utilities.f90 | 4 +- src/homogenization.f90 | 106 +++++++++++++++----------------- src/mesh/FEM_utilities.f90 | 2 +- 4 files changed, 53 insertions(+), 62 deletions(-) diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 0f7590782..8f10d39b2 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -193,11 +193,10 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS else validCalculation if (debugCPFEM%extensive) print'(a,i8,1x,i2)', '<< CPFEM >> calculation for elFE ip ',elFE,ip - call homogenization_mechanical_response(dt,[ip,ip],[elCP,elCP]) + call homogenization_mechanical_response(dt,(elCP-1)*discretization_nIPs + ip,(elCP-1)*discretization_nIPs + ip) if (.not. terminallyIll) & call homogenization_mechanical_response2(dt,[ip,ip],[elCP,elCP]) - terminalIllness: if (terminallyIll) then call random_number(rnd) diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index d95580145..34a976644 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -812,9 +812,9 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& homogenization_F = reshape(F,[3,3,product(cells(1:2))*cells3]) ! set materialpoint target F to estimated field - call homogenization_mechanical_response(Delta_t,[1,1],[1,product(cells(1:2))*cells3]) ! calculate P field + call homogenization_mechanical_response(Delta_t,1,product(cells(1:2))*cells3) ! calculate P field if (.not. terminallyIll) & - call homogenization_thermal_response(Delta_t,[1,1],[1,product(cells(1:2))*cells3]) + call homogenization_thermal_response(Delta_t,1,product(cells(1:2))*cells3) if (.not. terminallyIll) & call homogenization_mechanical_response2(Delta_t,[1,1],[1,product(cells(1:2))*cells3]) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 8c9ecaecb..8539df994 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -222,14 +222,13 @@ end subroutine homogenization_init !-------------------------------------------------------------------------------------------------- !> @brief !-------------------------------------------------------------------------------------------------- -subroutine homogenization_mechanical_response(Delta_t,FEsolving_execIP,FEsolving_execElem) +subroutine homogenization_mechanical_response(Delta_t,cell_start,cell_end) real(pReal), intent(in) :: Delta_t !< time increment - integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP + integer, intent(in) :: & + cell_start, cell_end integer :: & NiterationMPstate, & - ip, & !< integration point number - el, & !< element number co, ce, ho, en logical :: & converged @@ -237,45 +236,42 @@ subroutine homogenization_mechanical_response(Delta_t,FEsolving_execIP,FEsolving doneAndHappy - !$OMP PARALLEL DO PRIVATE(ce,co,en,ho,NiterationMPstate,converged,doneAndHappy) - do el = FEsolving_execElem(1),FEsolving_execElem(2) + !$OMP PARALLEL DO PRIVATE(en,ho,co,NiterationMPstate,converged,doneAndHappy) + do ce = cell_start, cell_end - do ip = FEsolving_execIP(1),FEsolving_execIP(2) - ce = (el-1)*discretization_nIPs + ip - en = material_homogenizationEntry(ce) - ho = material_homogenizationID(ce) + en = material_homogenizationEntry(ce) + ho = material_homogenizationID(ce) - call phase_restore(ce,.false.) ! wrong name (is more a forward function) + call phase_restore(ce,.false.) ! wrong name (is more a forward function) - if(homogState(ho)%sizeState > 0) homogState(ho)%state(:,en) = homogState(ho)%state0(:,en) - if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%state(:,en) = damageState_h(ho)%state0(:,en) - call damage_partition(ce) + if(homogState(ho)%sizeState > 0) homogState(ho)%state(:,en) = homogState(ho)%state0(:,en) + if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%state(:,en) = damageState_h(ho)%state0(:,en) + call damage_partition(ce) - doneAndHappy = [.false.,.true.] + doneAndHappy = [.false.,.true.] - NiterationMPstate = 0 - convergenceLooping: do while (.not. (terminallyIll .or. doneAndHappy(1)) & - .and. NiterationMPstate < num%nMPstate) - NiterationMPstate = NiterationMPstate + 1 + NiterationMPstate = 0 + convergenceLooping: do while (.not. (terminallyIll .or. doneAndHappy(1)) & + .and. NiterationMPstate < num%nMPstate) + NiterationMPstate = NiterationMPstate + 1 - call mechanical_partition(homogenization_F(1:3,1:3,ce),ce) - converged = all([(phase_mechanical_constitutive(Delta_t,co,ce),co=1,homogenization_Nconstituents(ho))]) - if (converged) then - doneAndHappy = mechanical_updateState(Delta_t,homogenization_F(1:3,1:3,ce),ce) - converged = all(doneAndHappy) - else - doneAndHappy = [.true.,.false.] - endif - enddo convergenceLooping + call mechanical_partition(homogenization_F(1:3,1:3,ce),ce) + converged = all([(phase_mechanical_constitutive(Delta_t,co,ce),co=1,homogenization_Nconstituents(ho))]) + if (converged) then + doneAndHappy = mechanical_updateState(Delta_t,homogenization_F(1:3,1:3,ce),ce) + converged = all(doneAndHappy) + else + doneAndHappy = [.true.,.false.] + end if + end do convergenceLooping - converged = converged .and. all([(phase_damage_constitutive(Delta_t,co,ce),co=1,homogenization_Nconstituents(ho))]) + converged = converged .and. all([(phase_damage_constitutive(Delta_t,co,ce),co=1,homogenization_Nconstituents(ho))]) - if (.not. converged) then - if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill' - terminallyIll = .true. - endif - enddo - enddo + if (.not. converged) then + if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill' + terminallyIll = .true. + end if + end do !$OMP END PARALLEL DO end subroutine homogenization_mechanical_response @@ -284,31 +280,27 @@ end subroutine homogenization_mechanical_response !-------------------------------------------------------------------------------------------------- !> @brief !-------------------------------------------------------------------------------------------------- -subroutine homogenization_thermal_response(Delta_t,FEsolving_execIP,FEsolving_execElem) +subroutine homogenization_thermal_response(Delta_t,cell_start,cell_end) real(pReal), intent(in) :: Delta_t !< time increment - integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP + integer, intent(in) :: & + cell_start, cell_end integer :: & - ip, & !< integration point number - el, & !< element number co, ce, ho - !$OMP PARALLEL DO PRIVATE(ho,ce) - do el = FEsolving_execElem(1),FEsolving_execElem(2) + !$OMP PARALLEL DO PRIVATE(ho) + do ce = cell_start, cell_end if (terminallyIll) continue - do ip = FEsolving_execIP(1),FEsolving_execIP(2) - ce = (el-1)*discretization_nIPs + ip - ho = material_homogenizationID(ce) - call thermal_partition(ce) - do co = 1, homogenization_Nconstituents(ho) - if (.not. phase_thermal_constitutive(Delta_t,material_phaseID(co,ce),material_phaseEntry(co,ce))) then - if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill' - terminallyIll = .true. - endif - enddo - enddo - enddo + ho = material_homogenizationID(ce) + call thermal_partition(ce) + do co = 1, homogenization_Nconstituents(ho) + if (.not. phase_thermal_constitutive(Delta_t,material_phaseID(co,ce),material_phaseEntry(co,ce))) then + if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill' + terminallyIll = .true. + end if + end do + end do !$OMP END PARALLEL DO end subroutine homogenization_thermal_response @@ -331,11 +323,11 @@ subroutine homogenization_mechanical_response2(Delta_t,FEsolving_execIP,FEsolvin elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2) IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2) ce = (el-1)*discretization_nIPs + ip - ho = material_homogenizationID(ce) - do co = 1, homogenization_Nconstituents(ho) - call crystallite_orientations(co,ip,el) + ho = material_homogenizationID(ce) + do co = 1, homogenization_Nconstituents(ho) + call crystallite_orientations(co,ip,el) enddo - call mechanical_homogenize(Delta_t,ce) + call mechanical_homogenize(Delta_t,ce) enddo IpLooping3 enddo elementLooping3 !$OMP END PARALLEL DO diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 703c6c818..4e31a0475 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -150,7 +150,7 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData) print'(/,1x,a)', '... evaluating constitutive response ......................................' - call homogenization_mechanical_response(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) ! calculate P field + call homogenization_mechanical_response(timeinc,1,mesh_maxNips*mesh_NcpElems) ! calculate P field if (.not. terminallyIll) & call homogenization_mechanical_response2(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) cutBack = .false. From 5f0a630fa6000a236bb359c83603cd74966b6d52 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 5 Feb 2022 08:12:35 +0100 Subject: [PATCH 11/13] remove deprecated phaseAt/phaseMemberAt --- src/phase_mechanical_plastic_nonlocal.f90 | 26 +++++++++++------------ 1 file changed, 12 insertions(+), 14 deletions(-) diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index dec7ea8e7..9d809e0dc 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -659,10 +659,10 @@ module subroutine nonlocal_dependentState(ph, en) do n = 1,nIPneighbors neighbor_el = geom(ph)%IPneighborhood(1,n,en) neighbor_ip = geom(ph)%IPneighborhood(2,n,en) - no = material_phasememberAt(1,neighbor_ip,neighbor_el) - if (neighbor_el > 0 .and. neighbor_ip > 0) then - if (material_phaseAt(1,neighbor_el) == ph) then + if (neighbor_el > 0 .and. neighbor_ip > 0) then + if (material_phaseID(1,(neighbor_el-1)*discretization_nIPs + neighbor_ip) == ph) then + no = material_phaseEntry(1,(neighbor_el-1)*discretization_nIPs + neighbor_ip) nRealNeighbors = nRealNeighbors + 1.0_pReal rho_neighbor0 = getRho0(ph,no) @@ -1145,7 +1145,6 @@ function rhoDotFlux(timestep,ph,en) en integer :: & - neighbor_ph, & !< phase of my neighbor's plasticity ns, & !< short notation for the total number of active slip systems c, & !< character of dislocation n, & !< index of my current neighbor @@ -1251,8 +1250,8 @@ function rhoDotFlux(timestep,ph,en) neighbor_el = geom(ph)%IPneighborhood(1,n,en) neighbor_ip = geom(ph)%IPneighborhood(2,n,en) neighbor_n = geom(ph)%IPneighborhood(3,n,en) - np = material_phaseAt(1,neighbor_el) - no = material_phasememberAt(1,neighbor_ip,neighbor_el) + np = material_phaseID(1,(neighbor_el-1)*discretization_nIPs + neighbor_ip) + no = material_phaseEntry(1,(neighbor_el-1)*discretization_nIPs + neighbor_ip) opposite_neighbor = n + mod(n,2) - mod(n+1,2) opposite_el = geom(ph)%IPneighborhood(1,opposite_neighbor,en) @@ -1260,7 +1259,6 @@ function rhoDotFlux(timestep,ph,en) opposite_n = geom(ph)%IPneighborhood(3,opposite_neighbor,en) if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient - neighbor_ph = material_phaseAt(1,neighbor_el) neighbor_F = phase_mechanical_F(np)%data(1:3,1:3,no) neighbor_Fe = matmul(neighbor_F, math_inv33(phase_mechanical_Fp(np)%data(1:3,1:3,no))) Favg = 0.5_pReal * (my_F + neighbor_F) @@ -1279,12 +1277,12 @@ function rhoDotFlux(timestep,ph,en) !* The entering flux from my neighbor will be distributed on my slip systems according to the !* compatibility if (neighbor_n > 0) then - if (phase_plasticity(material_phaseAt(1,neighbor_el)) == PLASTIC_NONLOCAL_ID .and. & + if (phase_plasticity(np) == PLASTIC_NONLOCAL_ID .and. & any(dependentState(ph)%compatibility(:,:,:,n,en) > 0.0_pReal)) then forall (s = 1:ns, t = 1:4) - neighbor_v0(s,t) = plasticState(np)%state0(iV (s,t,neighbor_ph),no) - neighbor_rhoSgl0(s,t) = max(plasticState(np)%state0(iRhoU(s,t,neighbor_ph),no),0.0_pReal) + neighbor_v0(s,t) = plasticState(np)%state0(iV (s,t,np),no) + neighbor_rhoSgl0(s,t) = max(plasticState(np)%state0(iRhoU(s,t,np),no),0.0_pReal) endforall where (neighbor_rhoSgl0 * IPvolume(neighbor_ip,neighbor_el) ** 0.667_pReal < prm%rho_min & @@ -1325,7 +1323,7 @@ function rhoDotFlux(timestep,ph,en) !* In case of reduced transmissivity, part of the leaving flux is stored as dead dislocation density. !* That means for an interface of zero transmissivity the leaving flux is fully converted to dead dislocations. if (opposite_n > 0) then - if (phase_plasticity(material_phaseAt(1,opposite_el)) == PLASTIC_NONLOCAL_ID) then + if (phase_plasticity(np) == PLASTIC_NONLOCAL_ID) then normal_me2neighbor_defConf = math_det33(Favg) & * matmul(math_inv33(transpose(Favg)),geom(ph)%IPareaNormal(1:3,n,en)) ! normal of the interface in (average) deformed configuration (pointing en => neighbor) @@ -1400,7 +1398,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,ip,el) associate(prm => param(ph)) ns = prm%sum_N_sl - en = material_phaseMemberAt(1,ip,el) + en = material_phaseEntry(1,(el-1)*discretization_nIPs + ip) !*** start out fully compatible my_compatibility = 0.0_pReal forall(s1 = 1:ns) my_compatibility(:,s1,s1,:) = 1.0_pReal @@ -1408,8 +1406,8 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,ip,el) neighbors: do n = 1,nIPneighbors neighbor_e = IPneighborhood(1,n,ip,el) neighbor_i = IPneighborhood(2,n,ip,el) - neighbor_me = material_phaseMemberAt(1,neighbor_i,neighbor_e) - neighbor_phase = material_phaseAt(1,neighbor_e) + neighbor_me = material_phaseEntry(1,(neighbor_e-1)*discretization_nIPs + neighbor_i) + neighbor_phase = material_phaseID(1,(neighbor_e-1)*discretization_nIPs + neighbor_i) if (neighbor_e <= 0 .or. neighbor_i <= 0) then !* FREE SURFACE From 6c032e3ce6f01f93d15c7a643b17c117ca54bebf Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 5 Feb 2022 09:03:10 +0100 Subject: [PATCH 12/13] remove deprecated mappings almost done with having a consistent access pattern solver ====== grid: x,y,z; mesh: el,ip homogenization ============== interface to solver: ce internal: ho,en phase ===== interface to homogenization: co,ce internal: ph,en --- src/material.f90 | 81 ++++++++++++++++++++++-------------------------- src/phase.f90 | 4 +-- 2 files changed, 39 insertions(+), 46 deletions(-) diff --git a/src/material.f90 b/src/material.f90 index fd1531964..bf78a77a6 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -17,16 +17,17 @@ module material implicit none private - type :: tRotationContainer - type(tRotation), dimension(:), allocatable :: data - end type - type :: tTensorContainer + type, public :: tRotationContainer + type(tRotation), dimension(:), allocatable :: data + end type tRotationContainer + + type, public :: tTensorContainer real(pReal), dimension(:,:,:), allocatable :: data - end type + end type tTensorContainer - type(tRotationContainer), dimension(:), allocatable :: material_O_0 - type(tTensorContainer), dimension(:), allocatable :: material_F_i_0 + type(tRotationContainer), dimension(:), allocatable, public, protected :: material_O_0 + type(tTensorContainer), dimension(:), allocatable, public, protected :: material_F_i_0 integer, dimension(:), allocatable, public, protected :: & homogenization_Nconstituents !< number of grains in each homogenization @@ -37,20 +38,14 @@ module material material_name_phase, & !< name of each phase material_name_homogenization !< name of each homogenization - integer, dimension(:), allocatable, public, protected :: & ! (elem) - material_homogenizationID, & !< per cell TODO: material_ID_homogenization - material_homogenizationEntry !< per cell TODO: material_entry_homogenization - integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,elem) - material_phaseAt, & !< phase ID of each element TODO: remove - material_phaseID, & !< per (constituent,cell) TODO: material_ID_phase - material_phaseEntry !< per (constituent,cell) TODO: material_entry_phase - integer, dimension(:,:,:), allocatable, public, protected :: & ! (constituent,IP,elem) - material_phaseMemberAt !TODO: remove + integer, dimension(:), allocatable, public, protected :: & ! (cell) + material_homogenizationID, & ! TODO: rename to material_ID_homogenization + material_homogenizationEntry ! TODO: rename to material_entry_homogenization + integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,cell) + material_phaseID, & ! TODO: rename to material_ID_phase + material_phaseEntry ! TODO: rename to material_entry_phase + public :: & - tTensorContainer, & - tRotationContainer, & - material_F_i_0, & - material_O_0, & material_init contains @@ -97,11 +92,12 @@ subroutine parse() counterPhase, & counterHomogenization - real(pReal) :: & - frac + real(pReal) :: v integer :: & - el, ip, co, ma, & - h, ce + el, ip, & + ho, ph, & + co, ce, & + ma materials => config_material%get('material') phases => config_material%get('phase') @@ -118,51 +114,48 @@ subroutine parse() #endif allocate(homogenization_Nconstituents(homogenizations%length)) - do h=1, homogenizations%length - homogenization => homogenizations%get(h) - homogenization_Nconstituents(h) = homogenization%get_asInt('N_constituents') + do ho=1, homogenizations%length + homogenization => homogenizations%get(ho) + homogenization_Nconstituents(ho) = homogenization%get_asInt('N_constituents') end do homogenization_maxNconstituents = maxval(homogenization_Nconstituents) allocate(counterPhase(phases%length),source=0) allocate(counterHomogenization(homogenizations%length),source=0) - allocate(material_phaseAt(homogenization_maxNconstituents,discretization_Nelems),source=0) - allocate(material_phaseMemberAt(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems),source=0) - - allocate(material_homogenizationID(discretization_nIPs*discretization_Nelems),source=0) allocate(material_homogenizationEntry(discretization_nIPs*discretization_Nelems),source=0) + allocate(material_phaseID(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0) allocate(material_phaseEntry(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0) do el = 1, discretization_Nelems - material => materials%get(discretization_materialAt(el)) - constituents => material%get('constituents') + material => materials%get(discretization_materialAt(el)) + ho = homogenizations%getIndex(material%get_asString('homogenization')) do ip = 1, discretization_nIPs ce = (el-1)*discretization_nIPs + ip - material_homogenizationID(ce) = homogenizations%getIndex(material%get_asString('homogenization')) - counterHomogenization(material_homogenizationID(ce)) = counterHomogenization(material_homogenizationID(ce)) + 1 - material_homogenizationEntry(ce) = counterHomogenization(material_homogenizationID(ce)) + material_homogenizationID(ce) = ho + counterHomogenization(ho) = counterHomogenization(ho) + 1 + material_homogenizationEntry(ce) = counterHomogenization(ho) end do - frac = 0.0_pReal + v = 0.0_pReal + constituents => material%get('constituents') do co = 1, constituents%length constituent => constituents%get(co) - frac = frac + constituent%get_asFloat('v') + v = v + constituent%get_asFloat('v') - material_phaseAt(co,el) = phases%getIndex(constituent%get_asString('phase')) + ph = phases%getIndex(constituent%get_asString('phase')) do ip = 1, discretization_nIPs ce = (el-1)*discretization_nIPs + ip - counterPhase(material_phaseAt(co,el)) = counterPhase(material_phaseAt(co,el)) + 1 - material_phaseMemberAt(co,ip,el) = counterPhase(material_phaseAt(co,el)) - material_phaseEntry(co,ce) = counterPhase(material_phaseAt(co,el)) - material_phaseID(co,ce) = material_phaseAt(co,el) + material_phaseID(co,ce) = ph + counterPhase(ph) = counterPhase(ph) + 1 + material_phaseEntry(co,ce) = counterPhase(ph) end do end do - if (dNeq(frac,1.0_pReal,1.e-12_pReal)) call IO_error(153,ext_msg='constituent') + if (dNeq(v,1.0_pReal,1.e-12_pReal)) call IO_error(153,ext_msg='constituent') end do diff --git a/src/phase.f90 b/src/phase.f90 index 9943c777a..354c738f1 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -590,8 +590,8 @@ subroutine crystallite_orientations(co,ip,el) call phase_O(ph)%data(en)%fromMatrix(transpose(math_rotationalPart(mechanical_F_e(ph,en)))) - if (plasticState(material_phaseAt(1,el))%nonlocal) & - call plastic_nonlocal_updateCompatibility(phase_O,material_phaseAt(1,el),ip,el) + if (plasticState(material_phaseID(1,(el-1)*discretization_nIPs + ip))%nonlocal) & + call plastic_nonlocal_updateCompatibility(phase_O,material_phaseID(1,(el-1)*discretization_nIPs + ip),ip,el) end subroutine crystallite_orientations From 9386ba5ef5a453e338d4f1f5c189025bee042f64 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 6 Feb 2022 12:59:41 +0100 Subject: [PATCH 13/13] allocate not needed --- src/phase_mechanical.f90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index de9c2bea2..fe62e16f5 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -996,16 +996,15 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_) subLi0, & subF0, & subF - real(pReal), dimension(:), allocatable :: subState0 + real(pReal), dimension(plasticState(material_phaseID(co,ce))%sizeState) :: subState0 ph = material_phaseID(co,ce) en = material_phaseEntry(co,ce) - sizeDotState = plasticState(ph)%sizeDotState + subState0 = plasticState(ph)%state0(:,en) subLi0 = phase_mechanical_Li0(ph)%data(1:3,1:3,en) subLp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,en) - allocate(subState0,source=plasticState(ph)%State0(:,en)) subFp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,en) subFi0 = phase_mechanical_Fi0(ph)%data(1:3,1:3,en) subF0 = phase_mechanical_F0(ph)%data(1:3,1:3,en) @@ -1038,7 +1037,7 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_) subStep = num%subStepSizeCryst * subStep phase_mechanical_Fp(ph)%data(1:3,1:3,en) = subFp0 phase_mechanical_Fi(ph)%data(1:3,1:3,en) = subFi0 - phase_mechanical_S(ph)%data(1:3,1:3,en) = phase_mechanical_S0(ph)%data(1:3,1:3,en) ! why no subS0 ? is S0 of any use? + phase_mechanical_S(ph)%data(1:3,1:3,en) = phase_mechanical_S0(ph)%data(1:3,1:3,en) if (subStep < 1.0_pReal) then ! actual (not initial) cutback phase_mechanical_Lp(ph)%data(1:3,1:3,en) = subLp0 phase_mechanical_Li(ph)%data(1:3,1:3,en) = subLi0 @@ -1050,6 +1049,7 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_) !-------------------------------------------------------------------------------------------------- ! prepare for integration if (todo) then + sizeDotState = plasticState(ph)%sizeDotState subF = subF0 & + subStep * (phase_mechanical_F(ph)%data(1:3,1:3,en) - phase_mechanical_F0(ph)%data(1:3,1:3,en)) converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * Delta_t,ph,en)