From c46b4d90a613a2ad881d9eec578326fc21e918ed Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 29 Sep 2020 12:48:29 +0200 Subject: [PATCH 1/5] modularizing --- src/crystallite.f90 | 22 +++++++++++++++++----- 1 file changed, 17 insertions(+), 5 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 0bd816142..c0c140571 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1241,7 +1241,7 @@ subroutine integrateStateFPI(todo) enddo; enddo; enddo !$OMP END PARALLEL DO - if (nonlocalBroken) call nonlocalConvergenceCheck + call nonlocalConvergenceCheck contains @@ -1328,7 +1328,7 @@ subroutine integrateStateEuler(todo) enddo; enddo; enddo !$OMP END PARALLEL DO - if (nonlocalBroken) call nonlocalConvergenceCheck + call nonlocalConvergenceCheck end subroutine integrateStateEuler @@ -1422,7 +1422,7 @@ subroutine integrateStateAdaptiveEuler(todo) enddo; enddo; enddo !$OMP END PARALLEL DO - if (nonlocalBroken) call nonlocalConvergenceCheck + call nonlocalConvergenceCheck end subroutine integrateStateAdaptiveEuler @@ -1612,7 +1612,7 @@ subroutine integrateStateRK(todo,A,B,CC,DB) enddo; enddo; enddo !$OMP END PARALLEL DO - if(nonlocalBroken) call nonlocalConvergenceCheck + call nonlocalConvergenceCheck end subroutine integrateStateRK @@ -1624,7 +1624,19 @@ end subroutine integrateStateRK subroutine nonlocalConvergenceCheck integer :: e,i,p - + logical :: nonlocal_broken + + nonlocal_broken = .false. + !$OMP PARALLEL DO PRIVATE(p) + do e = FEsolving_execElem(1),FEsolving_execElem(2) + p = material_phaseAt(1,e) + do i = FEsolving_execIP(1),FEsolving_execIP(2) + if(plasticState(p)%nonlocal .and. .not. crystallite_converged(1,i,e)) nonlocal_broken = .true. + enddo + enddo + !$OMP END PARALLEL DO + + if(.not. nonlocal_broken) return !$OMP PARALLEL DO PRIVATE(p) do e = FEsolving_execElem(1),FEsolving_execElem(2) p = material_phaseAt(1,e) From f1e96489cc3d59fe11420a9d395caf025fb234f1 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 29 Sep 2020 12:56:12 +0200 Subject: [PATCH 2/5] better readable somehow on the cost of the nonlocal performance --- src/crystallite.f90 | 43 ++++++++++++++++--------------------------- 1 file changed, 16 insertions(+), 27 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index c0c140571..5f4b3bbf7 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1147,16 +1147,15 @@ subroutine integrateStateFPI(todo) plastic_dotState real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState logical :: & - nonlocalBroken, broken + broken + - nonlocalBroken = .false. !$OMP PARALLEL DO PRIVATE(size_pl,size_so,r,zeta,p,c,plastic_dotState,source_dotState,broken) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) - p = material_phaseAt(g,e) - if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then - + if(todo(g,i,e)) then + p = material_phaseAt(g,e) c = material_phaseMemberAt(g,i,e) broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & @@ -1164,7 +1163,6 @@ subroutine integrateStateFPI(todo) crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle size_pl = plasticState(p)%sizeDotState @@ -1236,7 +1234,6 @@ subroutine integrateStateFPI(todo) endif enddo iteration - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. endif enddo; enddo; enddo !$OMP END PARALLEL DO @@ -1284,16 +1281,15 @@ subroutine integrateStateEuler(todo) s, & sizeDotState logical :: & - nonlocalBroken, broken + broken - nonlocalBroken = .false. !$OMP PARALLEL DO PRIVATE (sizeDotState,p,c,broken) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) - p = material_phaseAt(g,e) - if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then + if(todo(g,i,e)) then + p = material_phaseAt(g,e) c = material_phaseMemberAt(g,i,e) broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & @@ -1301,7 +1297,6 @@ subroutine integrateStateEuler(todo) crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle sizeDotState = plasticState(p)%sizeDotState @@ -1318,11 +1313,9 @@ subroutine integrateStateEuler(todo) broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle broken = integrateStress(g,i,e) - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. crystallite_converged(g,i,e) = .not. broken endif enddo; enddo; enddo @@ -1349,20 +1342,19 @@ subroutine integrateStateAdaptiveEuler(todo) s, & sizeDotState logical :: & - nonlocalBroken, broken + broken real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: residuum_plastic real(pReal), dimension(constitutive_source_maxSizeDotState,maxval(phase_Nsources)) :: residuum_source - nonlocalBroken = .false. !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,residuum_plastic,residuum_source,broken) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) broken = .false. - p = material_phaseAt(g,e) - if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then + if(todo(g,i,e)) then + p = material_phaseAt(g,e) c = material_phaseMemberAt(g,i,e) broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & @@ -1418,7 +1410,6 @@ subroutine integrateStateAdaptiveEuler(todo) enddo endif - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. enddo; enddo; enddo !$OMP END PARALLEL DO @@ -1503,19 +1494,18 @@ subroutine integrateStateRK(todo,A,B,CC,DB) s, & sizeDotState logical :: & - nonlocalBroken, broken + broken real(pReal), dimension(constitutive_source_maxSizeDotState,size(B),maxval(phase_Nsources)) :: source_RKdotState real(pReal), dimension(constitutive_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState - nonlocalBroken = .false. !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,plastic_RKdotState,source_RKdotState,broken) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) broken = .false. - p = material_phaseAt(g,e) - if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then + if(todo(g,i,e)) then + p = material_phaseAt(g,e) c = material_phaseMemberAt(g,i,e) broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & @@ -1608,7 +1598,6 @@ subroutine integrateStateRK(todo,A,B,CC,DB) crystallite_converged(g,i,e) = .not. broken endif - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. enddo; enddo; enddo !$OMP END PARALLEL DO @@ -1624,8 +1613,8 @@ end subroutine integrateStateRK subroutine nonlocalConvergenceCheck integer :: e,i,p - logical :: nonlocal_broken - + logical :: nonlocal_broken + nonlocal_broken = .false. !$OMP PARALLEL DO PRIVATE(p) do e = FEsolving_execElem(1),FEsolving_execElem(2) @@ -1635,7 +1624,7 @@ subroutine nonlocalConvergenceCheck enddo enddo !$OMP END PARALLEL DO - + if(.not. nonlocal_broken) return !$OMP PARALLEL DO PRIVATE(p) do e = FEsolving_execElem(1),FEsolving_execElem(2) From 587d5ee445394808eb55ad2e9164575428e589ea Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 29 Sep 2020 13:13:53 +0200 Subject: [PATCH 3/5] no need for two loops --- src/crystallite.f90 | 97 ++++++++++++--------------------------------- 1 file changed, 25 insertions(+), 72 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 5f4b3bbf7..d21fd453f 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -460,16 +460,18 @@ function crystallite_stress() math_inv33(crystallite_Fi(1:3,1:3,c,i,e))) crystallite_subdt(c,i,e) = crystallite_subStep(c,i,e) * crystallite_dt(c,i,e) crystallite_converged(c,i,e) = .false. + call integrateState(c,i,e) endif enddo enddo enddo elementLooping3 !$OMP END PARALLEL DO + + call nonlocalConvergenceCheck !-------------------------------------------------------------------------------------------------- ! integrate --- requires fully defined state array (basic + dependent state) - if (any(todo)) call integrateState(todo) ! TODO: unroll into proper elementloop to avoid N^2 for single point evaluation where(.not. crystallite_converged .and. crystallite_subStep > num%subStepMinCryst) & ! do not try non-converged but fully cutbacked any further todo = .true. ! TODO: again unroll this into proper elementloop to avoid N^2 for single point evaluation @@ -1125,9 +1127,8 @@ end function integrateStress !> @brief integrate stress, state with adaptive 1st order explicit Euler method !> using Fixed Point Iteration to adapt the stepsize !-------------------------------------------------------------------------------------------------- -subroutine integrateStateFPI(todo) +subroutine integrateStateFPI(g,i,e) - logical, dimension(:,:,:), intent(in) :: todo integer :: & NiterationState, & !< number of iterations in state loop e, & !< element index in element loop @@ -1150,11 +1151,6 @@ subroutine integrateStateFPI(todo) broken - !$OMP PARALLEL DO PRIVATE(size_pl,size_so,r,zeta,p,c,plastic_dotState,source_dotState,broken) - do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1),FEsolving_execIP(2) - do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) - if(todo(g,i,e)) then p = material_phaseAt(g,e) c = material_phaseMemberAt(g,i,e) @@ -1163,7 +1159,7 @@ subroutine integrateStateFPI(todo) crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken) cycle + if(broken) return size_pl = plasticState(p)%sizeDotState plasticState(p)%state(1:size_pl,c) = plasticState(p)%subState0(1:size_pl,c) & @@ -1234,11 +1230,7 @@ subroutine integrateStateFPI(todo) endif enddo iteration - endif - enddo; enddo; enddo - !$OMP END PARALLEL DO - call nonlocalConvergenceCheck contains @@ -1268,9 +1260,7 @@ end subroutine integrateStateFPI !-------------------------------------------------------------------------------------------------- !> @brief integrate state with 1st order explicit Euler method !-------------------------------------------------------------------------------------------------- -subroutine integrateStateEuler(todo) - - logical, dimension(:,:,:), intent(in) :: todo +subroutine integrateStateEuler(g,i,e) integer :: & e, & !< element index in element loop @@ -1283,11 +1273,6 @@ subroutine integrateStateEuler(todo) logical :: & broken - !$OMP PARALLEL DO PRIVATE (sizeDotState,p,c,broken) - do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1),FEsolving_execIP(2) - do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) - if(todo(g,i,e)) then p = material_phaseAt(g,e) c = material_phaseMemberAt(g,i,e) @@ -1297,7 +1282,7 @@ subroutine integrateStateEuler(todo) crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken) cycle + if(broken) return sizeDotState = plasticState(p)%sizeDotState plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & @@ -1313,15 +1298,10 @@ subroutine integrateStateEuler(todo) broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) - if(broken) cycle + if(broken) return broken = integrateStress(g,i,e) crystallite_converged(g,i,e) = .not. broken - endif - enddo; enddo; enddo - !$OMP END PARALLEL DO - - call nonlocalConvergenceCheck end subroutine integrateStateEuler @@ -1329,9 +1309,7 @@ end subroutine integrateStateEuler !-------------------------------------------------------------------------------------------------- !> @brief integrate stress, state with 1st order Euler method with adaptive step size !-------------------------------------------------------------------------------------------------- -subroutine integrateStateAdaptiveEuler(todo) - - logical, dimension(:,:,:), intent(in) :: todo +subroutine integrateStateAdaptiveEuler(g,i,e) integer :: & e, & ! element index in element loop @@ -1347,13 +1325,7 @@ subroutine integrateStateAdaptiveEuler(todo) real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: residuum_plastic real(pReal), dimension(constitutive_source_maxSizeDotState,maxval(phase_Nsources)) :: residuum_source - !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,residuum_plastic,residuum_source,broken) - do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1),FEsolving_execIP(2) - do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) - broken = .false. - if(todo(g,i,e)) then p = material_phaseAt(g,e) c = material_phaseMemberAt(g,i,e) @@ -1362,7 +1334,7 @@ subroutine integrateStateAdaptiveEuler(todo) crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken) cycle + if(broken) return sizeDotState = plasticState(p)%sizeDotState @@ -1381,17 +1353,17 @@ subroutine integrateStateAdaptiveEuler(todo) broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) - if(broken) cycle + if(broken) return broken = integrateStress(g,i,e) - if(broken) cycle + if(broken) return broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken) cycle + if(broken) return sizeDotState = plasticState(p)%sizeDotState @@ -1407,13 +1379,7 @@ subroutine integrateStateAdaptiveEuler(todo) + 0.5_pReal*sourceState(p)%p(s)%dotState(:,c)*crystallite_subdt(g,i,e), & sourceState(p)%p(s)%state(1:sizeDotState,c), & sourceState(p)%p(s)%atol(1:sizeDotState)) - enddo - - endif - enddo; enddo; enddo - !$OMP END PARALLEL DO - - call nonlocalConvergenceCheck + enddo end subroutine integrateStateAdaptiveEuler @@ -1421,9 +1387,9 @@ end subroutine integrateStateAdaptiveEuler !--------------------------------------------------------------------------------------------------- !> @brief Integrate state (including stress integration) with the classic Runge Kutta method !--------------------------------------------------------------------------------------------------- -subroutine integrateStateRK4(todo) +subroutine integrateStateRK4(g,i,e) - logical, dimension(:,:,:), intent(in) :: todo + integer :: g,i,e real(pReal), dimension(3,3), parameter :: & A = reshape([& @@ -1436,7 +1402,7 @@ subroutine integrateStateRK4(todo) real(pReal), dimension(4), parameter :: & 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] - call integrateStateRK(todo,A,B,C) + call integrateStateRK(g,i,e,A,B,C) end subroutine integrateStateRK4 @@ -1444,9 +1410,9 @@ end subroutine integrateStateRK4 !--------------------------------------------------------------------------------------------------- !> @brief Integrate state (including stress integration) with the Cash-Carp method !--------------------------------------------------------------------------------------------------- -subroutine integrateStateRKCK45(todo) +subroutine integrateStateRKCK45(g,i,e) - logical, dimension(:,:,:), intent(in) :: todo + integer :: g,i,e real(pReal), dimension(5,5), parameter :: & A = reshape([& @@ -1466,7 +1432,7 @@ subroutine integrateStateRKCK45(todo) [2825.0_pReal/27648.0_pReal, .0_pReal, 18575.0_pReal/48384.0_pReal,& 13525.0_pReal/55296.0_pReal, 277.0_pReal/14336.0_pReal, 1._pReal/4._pReal] - call integrateStateRK(todo,A,B,C,DB) + call integrateStateRK(g,i,e,A,B,C,DB) end subroutine integrateStateRKCK45 @@ -1475,9 +1441,8 @@ end subroutine integrateStateRKCK45 !> @brief Integrate state (including stress integration) with an explicit Runge-Kutta method or an !! embedded explicit Runge-Kutta method !-------------------------------------------------------------------------------------------------- -subroutine integrateStateRK(todo,A,B,CC,DB) +subroutine integrateStateRK(g,i,e,A,B,CC,DB) - logical, dimension(:,:,:), intent(in) :: todo real(pReal), dimension(:,:), intent(in) :: A real(pReal), dimension(:), intent(in) :: B, CC @@ -1498,13 +1463,6 @@ subroutine integrateStateRK(todo,A,B,CC,DB) real(pReal), dimension(constitutive_source_maxSizeDotState,size(B),maxval(phase_Nsources)) :: source_RKdotState real(pReal), dimension(constitutive_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState - !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,plastic_RKdotState,source_RKdotState,broken) - do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1),FEsolving_execIP(2) - do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) - broken = .false. - - if(todo(g,i,e)) then p = material_phaseAt(g,e) c = material_phaseMemberAt(g,i,e) @@ -1513,7 +1471,7 @@ subroutine integrateStateRK(todo,A,B,CC,DB) crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken) cycle + if(broken) return do stage = 1,size(A,1) sizeDotState = plasticState(p)%sizeDotState @@ -1558,7 +1516,7 @@ subroutine integrateStateRK(todo,A,B,CC,DB) if(broken) exit enddo - if(broken) cycle + if(broken) return sizeDotState = plasticState(p)%sizeDotState @@ -1587,21 +1545,16 @@ subroutine integrateStateRK(todo,A,B,CC,DB) sourceState(p)%p(s)%state(1:sizeDotState,c), & sourceState(p)%p(s)%atol(1:sizeDotState)) enddo - if(broken) cycle + if(broken) return broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) - if(broken) cycle + if(broken) return broken = integrateStress(g,i,e) crystallite_converged(g,i,e) = .not. broken - endif - enddo; enddo; enddo - !$OMP END PARALLEL DO - - call nonlocalConvergenceCheck end subroutine integrateStateRK From d0df748fc142481f29cdda9232634aade01b56f8 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 29 Sep 2020 13:36:15 +0200 Subject: [PATCH 4/5] cleaning --- src/crystallite.f90 | 474 ++++++++++++++++++++++---------------------- 1 file changed, 238 insertions(+), 236 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index d21fd453f..9411cdfc7 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -467,7 +467,7 @@ function crystallite_stress() enddo enddo elementLooping3 !$OMP END PARALLEL DO - + call nonlocalConvergenceCheck !-------------------------------------------------------------------------------------------------- @@ -1129,11 +1129,12 @@ end function integrateStress !-------------------------------------------------------------------------------------------------- subroutine integrateStateFPI(g,i,e) - integer :: & - NiterationState, & !< number of iterations in state loop + integer, intent(in) :: & e, & !< element index in element loop i, & !< integration point index in ip loop - g, & !< grain index in grain loop + g !< grain index in grain loop + integer :: & + NiterationState, & !< number of iterations in state loop p, & c, & s, & @@ -1150,86 +1151,85 @@ subroutine integrateStateFPI(g,i,e) logical :: & broken + p = material_phaseAt(g,e) + c = material_phaseMemberAt(g,i,e) - p = material_phaseAt(g,e) - c = material_phaseMemberAt(g,i,e) + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partionedF0, & + crystallite_Fi(1:3,1:3,g,i,e), & + crystallite_partionedFp0, & + crystallite_subdt(g,i,e), g,i,e,p,c) + if(broken) return - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partionedF0, & - crystallite_Fi(1:3,1:3,g,i,e), & - crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken) return + size_pl = plasticState(p)%sizeDotState + plasticState(p)%state(1:size_pl,c) = plasticState(p)%subState0(1:size_pl,c) & + + plasticState(p)%dotState (1:size_pl,c) & + * crystallite_subdt(g,i,e) + plastic_dotState(1:size_pl,2) = 0.0_pReal + do s = 1, phase_Nsources(p) + size_so(s) = sourceState(p)%p(s)%sizeDotState + sourceState(p)%p(s)%state(1:size_so(s),c) = sourceState(p)%p(s)%subState0(1:size_so(s),c) & + + sourceState(p)%p(s)%dotState (1:size_so(s),c) & + * crystallite_subdt(g,i,e) + source_dotState(1:size_so(s),2,s) = 0.0_pReal + enddo - size_pl = plasticState(p)%sizeDotState - plasticState(p)%state(1:size_pl,c) = plasticState(p)%subState0(1:size_pl,c) & - + plasticState(p)%dotState (1:size_pl,c) & - * crystallite_subdt(g,i,e) - plastic_dotState(1:size_pl,2) = 0.0_pReal - do s = 1, phase_Nsources(p) - size_so(s) = sourceState(p)%p(s)%sizeDotState - sourceState(p)%p(s)%state(1:size_so(s),c) = sourceState(p)%p(s)%subState0(1:size_so(s),c) & - + sourceState(p)%p(s)%dotState (1:size_so(s),c) & - * crystallite_subdt(g,i,e) - source_dotState(1:size_so(s),2,s) = 0.0_pReal - enddo + iteration: do NiterationState = 1, num%nState - iteration: do NiterationState = 1, num%nState + if(nIterationState > 1) plastic_dotState(1:size_pl,2) = plastic_dotState(1:size_pl,1) + plastic_dotState(1:size_pl,1) = plasticState(p)%dotState(:,c) + do s = 1, phase_Nsources(p) + if(nIterationState > 1) source_dotState(1:size_so(s),2,s) = source_dotState(1:size_so(s),1,s) + source_dotState(1:size_so(s),1,s) = sourceState(p)%p(s)%dotState(:,c) + enddo - if(nIterationState > 1) plastic_dotState(1:size_pl,2) = plastic_dotState(1:size_pl,1) - plastic_dotState(1:size_pl,1) = plasticState(p)%dotState(:,c) - do s = 1, phase_Nsources(p) - if(nIterationState > 1) source_dotState(1:size_so(s),2,s) = source_dotState(1:size_so(s),1,s) - source_dotState(1:size_so(s),1,s) = sourceState(p)%p(s)%dotState(:,c) - enddo + broken = integrateStress(g,i,e) + if(broken) exit iteration - broken = integrateStress(g,i,e) - if(broken) exit iteration + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partionedF0, & + crystallite_Fi(1:3,1:3,g,i,e), & + crystallite_partionedFp0, & + crystallite_subdt(g,i,e), g,i,e,p,c) + if(broken) exit iteration - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partionedF0, & - crystallite_Fi(1:3,1:3,g,i,e), & - crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken) exit iteration + zeta = damper(plasticState(p)%dotState(:,c),plastic_dotState(1:size_pl,1),& + plastic_dotState(1:size_pl,2)) + plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * zeta & + + plastic_dotState(1:size_pl,1) * (1.0_pReal - zeta) + r(1:size_pl) = plasticState(p)%state (1:size_pl,c) & + - plasticState(p)%subState0(1:size_pl,c) & + - plasticState(p)%dotState (1:size_pl,c) * crystallite_subdt(g,i,e) + plasticState(p)%state(1:size_pl,c) = plasticState(p)%state(1:size_pl,c) & + - r(1:size_pl) + crystallite_converged(g,i,e) = converged(r(1:size_pl), & + plasticState(p)%state(1:size_pl,c), & + plasticState(p)%atol(1:size_pl)) + do s = 1, phase_Nsources(p) + zeta = damper(sourceState(p)%p(s)%dotState(:,c), & + source_dotState(1:size_so(s),1,s),& + source_dotState(1:size_so(s),2,s)) + sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) * zeta & + + source_dotState(1:size_so(s),1,s)* (1.0_pReal - zeta) + r(1:size_so(s)) = sourceState(p)%p(s)%state (1:size_so(s),c) & + - sourceState(p)%p(s)%subState0(1:size_so(s),c) & + - sourceState(p)%p(s)%dotState (1:size_so(s),c) * crystallite_subdt(g,i,e) + sourceState(p)%p(s)%state(1:size_so(s),c) = sourceState(p)%p(s)%state(1:size_so(s),c) & + - r(1:size_so(s)) + crystallite_converged(g,i,e) = & + crystallite_converged(g,i,e) .and. converged(r(1:size_so(s)), & + sourceState(p)%p(s)%state(1:size_so(s),c), & + sourceState(p)%p(s)%atol(1:size_so(s))) + enddo - zeta = damper(plasticState(p)%dotState(:,c),plastic_dotState(1:size_pl,1),& - plastic_dotState(1:size_pl,2)) - plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * zeta & - + plastic_dotState(1:size_pl,1) * (1.0_pReal - zeta) - r(1:size_pl) = plasticState(p)%state (1:size_pl,c) & - - plasticState(p)%subState0(1:size_pl,c) & - - plasticState(p)%dotState (1:size_pl,c) * crystallite_subdt(g,i,e) - plasticState(p)%state(1:size_pl,c) = plasticState(p)%state(1:size_pl,c) & - - r(1:size_pl) - crystallite_converged(g,i,e) = converged(r(1:size_pl), & - plasticState(p)%state(1:size_pl,c), & - plasticState(p)%atol(1:size_pl)) - do s = 1, phase_Nsources(p) - zeta = damper(sourceState(p)%p(s)%dotState(:,c), & - source_dotState(1:size_so(s),1,s),& - source_dotState(1:size_so(s),2,s)) - sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) * zeta & - + source_dotState(1:size_so(s),1,s)* (1.0_pReal - zeta) - r(1:size_so(s)) = sourceState(p)%p(s)%state (1:size_so(s),c) & - - sourceState(p)%p(s)%subState0(1:size_so(s),c) & - - sourceState(p)%p(s)%dotState (1:size_so(s),c) * crystallite_subdt(g,i,e) - sourceState(p)%p(s)%state(1:size_so(s),c) = sourceState(p)%p(s)%state(1:size_so(s),c) & - - r(1:size_so(s)) - crystallite_converged(g,i,e) = & - crystallite_converged(g,i,e) .and. converged(r(1:size_so(s)), & - sourceState(p)%p(s)%state(1:size_so(s),c), & - sourceState(p)%p(s)%atol(1:size_so(s))) - enddo + if(crystallite_converged(g,i,e)) then + broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_Fe(1:3,1:3,g,i,e), & + crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) + exit iteration + endif - if(crystallite_converged(g,i,e)) then - broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_Fe(1:3,1:3,g,i,e), & - crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) - exit iteration - endif - - enddo iteration + enddo iteration contains @@ -1262,10 +1262,11 @@ end subroutine integrateStateFPI !-------------------------------------------------------------------------------------------------- subroutine integrateStateEuler(g,i,e) - integer :: & + integer, intent(in) :: & e, & !< element index in element loop i, & !< integration point index in ip loop - g, & !< grain index in grain loop + g !< grain index in grain loop + integer :: & p, & c, & s, & @@ -1273,35 +1274,34 @@ subroutine integrateStateEuler(g,i,e) logical :: & broken + p = material_phaseAt(g,e) + c = material_phaseMemberAt(g,i,e) - p = material_phaseAt(g,e) - c = material_phaseMemberAt(g,i,e) + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partionedF0, & + crystallite_Fi(1:3,1:3,g,i,e), & + crystallite_partionedFp0, & + crystallite_subdt(g,i,e), g,i,e,p,c) + if(broken) return - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partionedF0, & - crystallite_Fi(1:3,1:3,g,i,e), & - crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken) return + sizeDotState = plasticState(p)%sizeDotState + plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & + + plasticState(p)%dotState (1:sizeDotState,c) & + * crystallite_subdt(g,i,e) + do s = 1, phase_Nsources(p) + sizeDotState = sourceState(p)%p(s)%sizeDotState + sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & + + sourceState(p)%p(s)%dotState (1:sizeDotState,c) & + * crystallite_subdt(g,i,e) + enddo - sizeDotState = plasticState(p)%sizeDotState - plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & - + plasticState(p)%dotState (1:sizeDotState,c) & - * crystallite_subdt(g,i,e) - do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState - sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & - + sourceState(p)%p(s)%dotState (1:sizeDotState,c) & - * crystallite_subdt(g,i,e) - enddo + broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_Fe(1:3,1:3,g,i,e), & + crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) + if(broken) return - broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_Fe(1:3,1:3,g,i,e), & - crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) - if(broken) return - - broken = integrateStress(g,i,e) - crystallite_converged(g,i,e) = .not. broken + broken = integrateStress(g,i,e) + crystallite_converged(g,i,e) = .not. broken end subroutine integrateStateEuler @@ -1311,10 +1311,11 @@ end subroutine integrateStateEuler !-------------------------------------------------------------------------------------------------- subroutine integrateStateAdaptiveEuler(g,i,e) + integer, intent(in) :: & + e, & !< element index in element loop + i, & !< integration point index in ip loop + g !< grain index in grain loop integer :: & - e, & ! element index in element loop - i, & ! integration point index in ip loop - g, & ! grain index in grain loop p, & c, & s, & @@ -1326,60 +1327,60 @@ subroutine integrateStateAdaptiveEuler(g,i,e) real(pReal), dimension(constitutive_source_maxSizeDotState,maxval(phase_Nsources)) :: residuum_source - p = material_phaseAt(g,e) - c = material_phaseMemberAt(g,i,e) + p = material_phaseAt(g,e) + c = material_phaseMemberAt(g,i,e) - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partionedF0, & - crystallite_Fi(1:3,1:3,g,i,e), & - crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken) return + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partionedF0, & + crystallite_Fi(1:3,1:3,g,i,e), & + crystallite_partionedFp0, & + crystallite_subdt(g,i,e), g,i,e,p,c) + if(broken) return - sizeDotState = plasticState(p)%sizeDotState + sizeDotState = plasticState(p)%sizeDotState - residuum_plastic(1:sizeDotState) = - plasticState(p)%dotstate(1:sizeDotState,c) * 0.5_pReal * crystallite_subdt(g,i,e) - plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & - + plasticState(p)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) - do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState + residuum_plastic(1:sizeDotState) = - plasticState(p)%dotstate(1:sizeDotState,c) * 0.5_pReal * crystallite_subdt(g,i,e) + plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & + + plasticState(p)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) + do s = 1, phase_Nsources(p) + sizeDotState = sourceState(p)%p(s)%sizeDotState - residuum_source(1:sizeDotState,s) = - sourceState(p)%p(s)%dotstate(1:sizeDotState,c) & - * 0.5_pReal * crystallite_subdt(g,i,e) - sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & - + sourceState(p)%p(s)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) - enddo + residuum_source(1:sizeDotState,s) = - sourceState(p)%p(s)%dotstate(1:sizeDotState,c) & + * 0.5_pReal * crystallite_subdt(g,i,e) + sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & + + sourceState(p)%p(s)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) + enddo - broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_Fe(1:3,1:3,g,i,e), & - crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) - if(broken) return + broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_Fe(1:3,1:3,g,i,e), & + crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) + if(broken) return - broken = integrateStress(g,i,e) - if(broken) return + broken = integrateStress(g,i,e) + if(broken) return - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partionedF0, & - crystallite_Fi(1:3,1:3,g,i,e), & - crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken) return + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partionedF0, & + crystallite_Fi(1:3,1:3,g,i,e), & + crystallite_partionedFp0, & + crystallite_subdt(g,i,e), g,i,e,p,c) + if(broken) return - sizeDotState = plasticState(p)%sizeDotState - crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState) & - + 0.5_pReal * plasticState(p)%dotState(:,c) * crystallite_subdt(g,i,e), & - plasticState(p)%state(1:sizeDotState,c), & - plasticState(p)%atol(1:sizeDotState)) + sizeDotState = plasticState(p)%sizeDotState + crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState) & + + 0.5_pReal * plasticState(p)%dotState(:,c) * crystallite_subdt(g,i,e), & + plasticState(p)%state(1:sizeDotState,c), & + plasticState(p)%atol(1:sizeDotState)) - do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState - crystallite_converged(g,i,e) = & - crystallite_converged(g,i,e) .and. converged(residuum_source(1:sizeDotState,s) & - + 0.5_pReal*sourceState(p)%p(s)%dotState(:,c)*crystallite_subdt(g,i,e), & - sourceState(p)%p(s)%state(1:sizeDotState,c), & - sourceState(p)%p(s)%atol(1:sizeDotState)) - enddo + do s = 1, phase_Nsources(p) + sizeDotState = sourceState(p)%p(s)%sizeDotState + crystallite_converged(g,i,e) = & + crystallite_converged(g,i,e) .and. converged(residuum_source(1:sizeDotState,s) & + + 0.5_pReal*sourceState(p)%p(s)%dotState(:,c)*crystallite_subdt(g,i,e), & + sourceState(p)%p(s)%state(1:sizeDotState,c), & + sourceState(p)%p(s)%atol(1:sizeDotState)) + enddo end subroutine integrateStateAdaptiveEuler @@ -1389,7 +1390,7 @@ end subroutine integrateStateAdaptiveEuler !--------------------------------------------------------------------------------------------------- subroutine integrateStateRK4(g,i,e) - integer :: g,i,e + integer, intent(in) :: g,i,e real(pReal), dimension(3,3), parameter :: & A = reshape([& @@ -1412,7 +1413,7 @@ end subroutine integrateStateRK4 !--------------------------------------------------------------------------------------------------- subroutine integrateStateRKCK45(g,i,e) - integer :: g,i,e + integer, intent(in) :: g,i,e real(pReal), dimension(5,5), parameter :: & A = reshape([& @@ -1448,10 +1449,11 @@ subroutine integrateStateRK(g,i,e,A,B,CC,DB) real(pReal), dimension(:), intent(in) :: B, CC real(pReal), dimension(:), intent(in), optional :: DB + integer, intent(in) :: & + e, & !< element index in element loop + i, & !< integration point index in ip loop + g !< grain index in grain loop integer :: & - e, & ! element index in element loop - i, & ! integration point index in ip loop - g, & ! grain index in grain loop stage, & ! stage index in integration stage loop n, & p, & @@ -1463,97 +1465,97 @@ subroutine integrateStateRK(g,i,e,A,B,CC,DB) real(pReal), dimension(constitutive_source_maxSizeDotState,size(B),maxval(phase_Nsources)) :: source_RKdotState real(pReal), dimension(constitutive_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState - p = material_phaseAt(g,e) - c = material_phaseMemberAt(g,i,e) + p = material_phaseAt(g,e) + c = material_phaseMemberAt(g,i,e) - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partionedF0, & - crystallite_Fi(1:3,1:3,g,i,e), & - crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken) return + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partionedF0, & + crystallite_Fi(1:3,1:3,g,i,e), & + crystallite_partionedFp0, & + crystallite_subdt(g,i,e), g,i,e,p,c) + if(broken) return - do stage = 1,size(A,1) - sizeDotState = plasticState(p)%sizeDotState - plastic_RKdotState(1:sizeDotState,stage) = plasticState(p)%dotState(:,c) - plasticState(p)%dotState(:,c) = A(1,stage) * plastic_RKdotState(1:sizeDotState,1) - do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState - source_RKdotState(1:sizeDotState,stage,s) = sourceState(p)%p(s)%dotState(:,c) - sourceState(p)%p(s)%dotState(:,c) = A(1,stage) * source_RKdotState(1:sizeDotState,1,s) - enddo + do stage = 1,size(A,1) + sizeDotState = plasticState(p)%sizeDotState + plastic_RKdotState(1:sizeDotState,stage) = plasticState(p)%dotState(:,c) + plasticState(p)%dotState(:,c) = A(1,stage) * plastic_RKdotState(1:sizeDotState,1) + do s = 1, phase_Nsources(p) + sizeDotState = sourceState(p)%p(s)%sizeDotState + source_RKdotState(1:sizeDotState,stage,s) = sourceState(p)%p(s)%dotState(:,c) + sourceState(p)%p(s)%dotState(:,c) = A(1,stage) * source_RKdotState(1:sizeDotState,1,s) + enddo - do n = 2, stage - sizeDotState = plasticState(p)%sizeDotState - plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) & - + A(n,stage) * plastic_RKdotState(1:sizeDotState,n) - do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState - sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) & - + A(n,stage) * source_RKdotState(1:sizeDotState,n,s) - enddo - enddo + do n = 2, stage + sizeDotState = plasticState(p)%sizeDotState + plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) & + + A(n,stage) * plastic_RKdotState(1:sizeDotState,n) + do s = 1, phase_Nsources(p) + sizeDotState = sourceState(p)%p(s)%sizeDotState + sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) & + + A(n,stage) * source_RKdotState(1:sizeDotState,n,s) + enddo + enddo - sizeDotState = plasticState(p)%sizeDotState - plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & - + plasticState(p)%dotState (1:sizeDotState,c) & - * crystallite_subdt(g,i,e) - do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState - sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & - + sourceState(p)%p(s)%dotState (1:sizeDotState,c) & - * crystallite_subdt(g,i,e) - enddo - - broken = integrateStress(g,i,e,CC(stage)) - if(broken) exit - - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partionedF0, & - crystallite_Fi(1:3,1:3,g,i,e), & - crystallite_partionedFp0, & - crystallite_subdt(g,i,e)*CC(stage), g,i,e,p,c) - if(broken) exit - - enddo - if(broken) return - - sizeDotState = plasticState(p)%sizeDotState - - plastic_RKdotState(1:sizeDotState,size(B)) = plasticState (p)%dotState(:,c) - plasticState(p)%dotState(:,c) = matmul(plastic_RKdotState(1:sizeDotState,1:size(B)),B) - plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & - + plasticState(p)%dotState (1:sizeDotState,c) & + sizeDotState = plasticState(p)%sizeDotState + plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & + + plasticState(p)%dotState (1:sizeDotState,c) & + * crystallite_subdt(g,i,e) + do s = 1, phase_Nsources(p) + sizeDotState = sourceState(p)%p(s)%sizeDotState + sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & + + sourceState(p)%p(s)%dotState (1:sizeDotState,c) & * crystallite_subdt(g,i,e) - if(present(DB)) & - broken = .not. converged( matmul(plastic_RKdotState(1:sizeDotState,1:size(DB)),DB) & - * crystallite_subdt(g,i,e), & - plasticState(p)%state(1:sizeDotState,c), & - plasticState(p)%atol(1:sizeDotState)) + enddo - do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState + broken = integrateStress(g,i,e,CC(stage)) + if(broken) exit - source_RKdotState(1:sizeDotState,size(B),s) = sourceState(p)%p(s)%dotState(:,c) - sourceState(p)%p(s)%dotState(:,c) = matmul(source_RKdotState(1:sizeDotState,1:size(B),s),B) - sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & - + sourceState(p)%p(s)%dotState (1:sizeDotState,c) & - * crystallite_subdt(g,i,e) - if(present(DB)) & - broken = broken .or. .not. converged(matmul(source_RKdotState(1:sizeDotState,1:size(DB),s),DB) & - * crystallite_subdt(g,i,e), & - sourceState(p)%p(s)%state(1:sizeDotState,c), & - sourceState(p)%p(s)%atol(1:sizeDotState)) - enddo - if(broken) return + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partionedF0, & + crystallite_Fi(1:3,1:3,g,i,e), & + crystallite_partionedFp0, & + crystallite_subdt(g,i,e)*CC(stage), g,i,e,p,c) + if(broken) exit - broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_Fe(1:3,1:3,g,i,e), & - crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) - if(broken) return + enddo + if(broken) return - broken = integrateStress(g,i,e) - crystallite_converged(g,i,e) = .not. broken + sizeDotState = plasticState(p)%sizeDotState + + plastic_RKdotState(1:sizeDotState,size(B)) = plasticState (p)%dotState(:,c) + plasticState(p)%dotState(:,c) = matmul(plastic_RKdotState(1:sizeDotState,1:size(B)),B) + plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & + + plasticState(p)%dotState (1:sizeDotState,c) & + * crystallite_subdt(g,i,e) + if(present(DB)) & + broken = .not. converged( matmul(plastic_RKdotState(1:sizeDotState,1:size(DB)),DB) & + * crystallite_subdt(g,i,e), & + plasticState(p)%state(1:sizeDotState,c), & + plasticState(p)%atol(1:sizeDotState)) + + do s = 1, phase_Nsources(p) + sizeDotState = sourceState(p)%p(s)%sizeDotState + + source_RKdotState(1:sizeDotState,size(B),s) = sourceState(p)%p(s)%dotState(:,c) + sourceState(p)%p(s)%dotState(:,c) = matmul(source_RKdotState(1:sizeDotState,1:size(B),s),B) + sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & + + sourceState(p)%p(s)%dotState (1:sizeDotState,c) & + * crystallite_subdt(g,i,e) + if(present(DB)) & + broken = broken .or. .not. converged(matmul(source_RKdotState(1:sizeDotState,1:size(DB),s),DB) & + * crystallite_subdt(g,i,e), & + sourceState(p)%p(s)%state(1:sizeDotState,c), & + sourceState(p)%p(s)%atol(1:sizeDotState)) + enddo + if(broken) return + + broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_Fe(1:3,1:3,g,i,e), & + crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) + if(broken) return + + broken = integrateStress(g,i,e) + crystallite_converged(g,i,e) = .not. broken end subroutine integrateStateRK From bd5d557fbb560d08ace9d3fa8b0d3c91832d7639 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 29 Sep 2020 23:32:17 +0200 Subject: [PATCH 5/5] untangling the spaghetti subLp and subLi are local variables --- src/crystallite.f90 | 50 +++++++++++++-------------------------------- 1 file changed, 14 insertions(+), 36 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 9411cdfc7..ad1c2d4a7 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -54,12 +54,10 @@ module crystallite ! crystallite_Lp0, & !< plastic velocitiy grad at start of FE inc crystallite_partionedLp0, & !< plastic velocity grad at start of homog inc - crystallite_subLp0,& !< plastic velocity grad at start of crystallite inc ! crystallite_Li, & !< current intermediate velocitiy grad (end of converged time step) crystallite_Li0, & !< intermediate velocitiy grad at start of FE inc crystallite_partionedLi0, & !< intermediate velocity grad at start of homog inc - crystallite_subLi0, & !< intermediate velocity grad at start of crystallite inc ! crystallite_S0, & !< 2nd Piola-Kirchhoff stress vector at start of FE inc crystallite_partionedS0 !< 2nd Piola-Kirchhoff stress vector at start of homog inc @@ -183,7 +181,6 @@ subroutine crystallite_init crystallite_Li,crystallite_Lp, & crystallite_subF,crystallite_subF0, & crystallite_subFp0,crystallite_subFi0, & - crystallite_subLi0,crystallite_subLp0, & source = crystallite_partionedF) allocate(crystallite_dPdF(3,3,3,3,cMax,iMax,eMax),source=0.0_pReal) @@ -326,34 +323,17 @@ function crystallite_stress() startIP, endIP, & s logical, dimension(homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: todo !ToDo: need to set some values to false for different Ngrains - todo = .false. + real(pReal), dimension(:,:,:,:,:), allocatable :: & + subLp0,& !< plastic velocity grad at start of crystallite inc + subLi0 !< intermediate velocity grad at start of crystallite inc + + + todo = .false. + + subLp0 = crystallite_partionedLp0 + subLi0 = crystallite_partionedLi0 + -#ifdef DEBUG - if (debugCrystallite%selective & - .and. FEsolving_execElem(1) <= debugCrystallite%element & - .and. debugCrystallite%element <= FEsolving_execElem(2)) then - print'(/,a,i8,1x,i2,1x,i3)', '<< CRYST stress >> boundary and initial values at el ip ipc ', & - debugCrystallite%element,debugCrystallite%ip, debugCrystallite%grain - print'(a,/,3(12x,3(f14.9,1x)/))', '<< CRYST stress >> F ', & - transpose(crystallite_partionedF(1:3,1:3,debugCrystallite%grain, & - debugCrystallite%ip,debugCrystallite%element)) - print'(a,/,3(12x,3(f14.9,1x)/))', '<< CRYST stress >> F0 ', & - transpose(crystallite_partionedF0(1:3,1:3,debugCrystallite%grain, & - debugCrystallite%ip,debugCrystallite%element)) - print'(a,/,3(12x,3(f14.9,1x)/))', '<< CRYST stress >> Fp0', & - transpose(crystallite_partionedFp0(1:3,1:3,debugCrystallite%grain, & - debugCrystallite%ip,debugCrystallite%element)) - print'(a,/,3(12x,3(f14.9,1x)/))', '<< CRYST stress >> Fi0', & - transpose(crystallite_partionedFi0(1:3,1:3,debugCrystallite%grain, & - debugCrystallite%ip,debugCrystallite%element)) - print'(a,/,3(12x,3(f14.9,1x)/))', '<< CRYST stress >> Lp0', & - transpose(crystallite_partionedLp0(1:3,1:3,debugCrystallite%grain, & - debugCrystallite%ip,debugCrystallite%element)) - print'(a,/,3(12x,3(f14.9,1x)/))', '<< CRYST stress >> Li0', & - transpose(crystallite_partionedLi0(1:3,1:3,debugCrystallite%grain, & - debugCrystallite%ip,debugCrystallite%element)) - endif -#endif !-------------------------------------------------------------------------------------------------- ! initialize to starting condition @@ -370,9 +350,7 @@ function crystallite_stress() sourceState(material_phaseAt(c,e))%p(s)%partionedState0(:,material_phaseMemberAt(c,i,e)) enddo crystallite_subFp0(1:3,1:3,c,i,e) = crystallite_partionedFp0(1:3,1:3,c,i,e) - crystallite_subLp0(1:3,1:3,c,i,e) = crystallite_partionedLp0(1:3,1:3,c,i,e) crystallite_subFi0(1:3,1:3,c,i,e) = crystallite_partionedFi0(1:3,1:3,c,i,e) - crystallite_subLi0(1:3,1:3,c,i,e) = crystallite_partionedLi0(1:3,1:3,c,i,e) crystallite_subF0(1:3,1:3,c,i,e) = crystallite_partionedF0(1:3,1:3,c,i,e) crystallite_subFrac(c,i,e) = 0.0_pReal crystallite_subStep(c,i,e) = 1.0_pReal/num%subStepSizeCryst @@ -415,8 +393,8 @@ function crystallite_stress() todo(c,i,e) = crystallite_subStep(c,i,e) > 0.0_pReal ! still time left to integrate on? if (todo(c,i,e)) then crystallite_subF0 (1:3,1:3,c,i,e) = crystallite_subF(1:3,1:3,c,i,e) - crystallite_subLp0(1:3,1:3,c,i,e) = crystallite_Lp (1:3,1:3,c,i,e) - crystallite_subLi0(1:3,1:3,c,i,e) = crystallite_Li (1:3,1:3,c,i,e) + subLp0(1:3,1:3,c,i,e) = crystallite_Lp (1:3,1:3,c,i,e) + subLi0(1:3,1:3,c,i,e) = crystallite_Li (1:3,1:3,c,i,e) crystallite_subFp0(1:3,1:3,c,i,e) = crystallite_Fp (1:3,1:3,c,i,e) crystallite_subFi0(1:3,1:3,c,i,e) = crystallite_Fi (1:3,1:3,c,i,e) plasticState( material_phaseAt(c,e))%subState0(:,material_phaseMemberAt(c,i,e)) & @@ -435,8 +413,8 @@ function crystallite_stress() crystallite_Fi (1:3,1:3,c,i,e) = crystallite_subFi0(1:3,1:3,c,i,e) crystallite_S (1:3,1:3,c,i,e) = crystallite_S0 (1:3,1:3,c,i,e) if (crystallite_subStep(c,i,e) < 1.0_pReal) then ! actual (not initial) cutback - crystallite_Lp (1:3,1:3,c,i,e) = crystallite_subLp0(1:3,1:3,c,i,e) - crystallite_Li (1:3,1:3,c,i,e) = crystallite_subLi0(1:3,1:3,c,i,e) + crystallite_Lp (1:3,1:3,c,i,e) = subLp0(1:3,1:3,c,i,e) + crystallite_Li (1:3,1:3,c,i,e) = subLi0(1:3,1:3,c,i,e) endif plasticState (material_phaseAt(c,e))%state( :,material_phaseMemberAt(c,i,e)) & = plasticState(material_phaseAt(c,e))%subState0(:,material_phaseMemberAt(c,i,e))