From 570c0421000461cb193f3032c31f86a01d78a7f5 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 07:18:37 +0200 Subject: [PATCH 01/23] update of dependent state always in conjunction with state integration --- src/crystallite.f90 | 35 +++-------------------------------- 1 file changed, 3 insertions(+), 32 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index fb12c439b..d657b5b55 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -849,6 +849,9 @@ logical function integrateStress(ipc,ip,el,timeFraction) F = crystallite_subF(1:3,1:3,ipc,ip,el) endif + call constitutive_dependentState(crystallite_partionedF(1:3,1:3,ipc,ip,el), & + crystallite_Fp(1:3,1:3,ipc,ip,el),ipc,ip,el) + Lpguess = crystallite_Lp(1:3,1:3,ipc,ip,el) ! take as first guess Liguess = crystallite_Li(1:3,1:3,ipc,ip,el) ! take as first guess @@ -1060,10 +1063,6 @@ subroutine integrateStateFPI source_dotState(1:sizeDotState,1,s) = sourceState(p)%p(s)%dotState(:,c) enddo - call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), & - crystallite_Fp(1:3,1:3,g,i,e), & - g, i, e) - crystallite_todo(g,i,e) = integrateStress(g,i,e) if(.not. crystallite_todo(g,i,e)) exit iteration @@ -1202,10 +1201,6 @@ subroutine integrateStateEuler nonlocalBroken = .true. if(.not. crystallite_todo(g,i,e)) cycle - call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), & - crystallite_Fp(1:3,1:3,g,i,e), & - g, i, e) - crystallite_todo(g,i,e) = integrateStress(g,i,e) if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & nonlocalBroken = .true. @@ -1281,10 +1276,6 @@ subroutine integrateStateAdaptiveEuler nonlocalBroken = .true. if(.not. crystallite_todo(g,i,e)) cycle - call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), & - crystallite_Fp(1:3,1:3,g,i,e), & - g, i, e) - crystallite_todo(g,i,e) = integrateStress(g,i,e) if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & nonlocalBroken = .true. @@ -1415,10 +1406,6 @@ subroutine integrateStateRK4 * crystallite_subdt(g,i,e) enddo - call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), & - crystallite_Fp(1:3,1:3,g,i,e), & - g, i, e) - crystallite_todo(g,i,e) = integrateStress(g,i,e,CC(stage)) if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & nonlocalBroken = .true. @@ -1466,14 +1453,6 @@ subroutine integrateStateRK4 nonlocalBroken = .true. if(.not. crystallite_todo(g,i,e)) cycle - call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), & - crystallite_Fp(1:3,1:3,g,i,e), & - g, i, e) - - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & - nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) cycle - crystallite_todo(g,i,e) = integrateStress(g,i,e) if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & nonlocalBroken = .true. @@ -1583,10 +1562,6 @@ subroutine integrateStateRKCK45 * crystallite_subdt(g,i,e) enddo - call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), & - crystallite_Fp(1:3,1:3,g,i,e), & - g, i, e) - crystallite_todo(g,i,e) = integrateStress(g,i,e,CC(stage)) if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & nonlocalBroken = .true. @@ -1644,10 +1619,6 @@ subroutine integrateStateRKCK45 nonlocalBroken = .true. if(.not. crystallite_todo(g,i,e)) cycle - call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), & - crystallite_Fp(1:3,1:3,g,i,e), & - g, i, e) - crystallite_todo(g,i,e) = integrateStress(g,i,e) if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & nonlocalBroken = .true. From 9e926f1545d2f7218d63b5a0472e5d2221141fbe Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 07:41:55 +0200 Subject: [PATCH 02/23] centralize test for error --- src/constitutive.f90 | 20 +++++++++----- src/crystallite.f90 | 62 ++++++++++++-------------------------------- 2 files changed, 30 insertions(+), 52 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 94e590ab3..4402820e8 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -709,7 +709,7 @@ end subroutine constitutive_hooke_SandItsTangents !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -subroutine constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el) +function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el) result(broken) integer, intent(in) :: & ipc, & !< component-ID of integration point @@ -727,19 +727,22 @@ subroutine constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, real(pReal), dimension(3,3) :: & Mp integer :: & + phase, & ho, & !< homogenization tme, & !< thermal member position i, & !< counter in source loop instance, of + logical :: broken ho = material_homogenizationAt(el) tme = thermalMapping(ho)%p(ip,el) of = material_phasememberAt(ipc,ip,el) - instance = phase_plasticityInstance(material_phaseAt(ipc,el)) + phase = material_phaseAt(ipc,el) + instance = phase_plasticityInstance(phase) Mp = matmul(matmul(transpose(Fi),Fi),S) - plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) + plasticityType: select case (phase_plasticity(phase)) case (PLASTICITY_ISOTROPIC_ID) plasticityType call plastic_isotropic_dotState (Mp,instance,of) @@ -760,10 +763,11 @@ subroutine constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, call plastic_nonlocal_dotState (Mp,FArray,FpArray,temperature(ho)%p(tme),subdt, & instance,of,ip,el) end select plasticityType + broken = any(IEEE_is_NaN(plasticState(phase)%dotState(:,of))) - SourceLoop: do i = 1, phase_Nsources(material_phaseAt(ipc,el)) + SourceLoop: do i = 1, phase_Nsources(phase) - sourceType: select case (phase_source(i,material_phaseAt(ipc,el))) + sourceType: select case (phase_source(i,phase)) case (SOURCE_damage_anisoBrittle_ID) sourceType call source_damage_anisoBrittle_dotState (S, ipc, ip, el) !< correct stress? @@ -775,13 +779,15 @@ subroutine constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, call source_damage_anisoDuctile_dotState ( ipc, ip, el) case (SOURCE_thermal_externalheat_ID) sourceType - call source_thermal_externalheat_dotState(material_phaseAt(ipc,el),of) + call source_thermal_externalheat_dotState(phase,of) end select sourceType + broken = broken .or. any(IEEE_is_NaN(sourceState(phase)%p(i)%dotState(:,of))) + enddo SourceLoop -end subroutine constitutive_collectDotState +end function constitutive_collectDotState !-------------------------------------------------------------------------------------------------- diff --git a/src/crystallite.f90 b/src/crystallite.f90 index d657b5b55..03de4dea3 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1027,15 +1027,12 @@ subroutine integrateStateFPI p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) - call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_todo(g,i,e) = .not. 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) - crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c))) - do s = 1, phase_Nsources(p) - crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c))) - enddo + if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & nonlocalBroken = .true. if(.not. crystallite_todo(g,i,e)) cycle @@ -1066,15 +1063,12 @@ subroutine integrateStateFPI crystallite_todo(g,i,e) = integrateStress(g,i,e) if(.not. crystallite_todo(g,i,e)) exit iteration - call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_todo(g,i,e) = .not. 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) - crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c))) - do s = 1, phase_Nsources(p) - crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c))) - enddo + if(.not. crystallite_todo(g,i,e)) exit iteration sizeDotState = plasticState(p)%sizeDotState @@ -1172,15 +1166,12 @@ subroutine integrateStateEuler p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) - call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_todo(g,i,e) = .not. 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) - crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c))) - do s = 1, phase_Nsources(p) - crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c))) - enddo + if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & nonlocalBroken = .true. if(.not. crystallite_todo(g,i,e)) cycle @@ -1244,15 +1235,12 @@ subroutine integrateStateAdaptiveEuler p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) - call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_todo(g,i,e) = .not. 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) - crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c))) - do s = 1, phase_Nsources(p) - crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c))) - enddo + if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & nonlocalBroken = .true. if(.not. crystallite_todo(g,i,e)) cycle @@ -1281,15 +1269,12 @@ subroutine integrateStateAdaptiveEuler nonlocalBroken = .true. if(.not. crystallite_todo(g,i,e)) cycle - call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_todo(g,i,e) = .not. 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) - crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c))) - do s = 1, phase_Nsources(p) - crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c))) - enddo + if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & nonlocalBroken = .true. if(.not. crystallite_todo(g,i,e)) cycle @@ -1361,15 +1346,12 @@ subroutine integrateStateRK4 p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) - call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_todo(g,i,e) = .not. 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) - crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c))) - do s = 1, phase_Nsources(p) - crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c))) - enddo + if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & nonlocalBroken = .true. if(.not. crystallite_todo(g,i,e)) cycle @@ -1411,15 +1393,12 @@ subroutine integrateStateRK4 nonlocalBroken = .true. if(.not. crystallite_todo(g,i,e)) exit - call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_todo(g,i,e) = .not. 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) - crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c))) - do s = 1, phase_Nsources(p) - crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c))) - enddo + if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & nonlocalBroken = .true. if(.not. crystallite_todo(g,i,e)) exit @@ -1517,15 +1496,12 @@ subroutine integrateStateRKCK45 p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) - call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_todo(g,i,e) = .not. 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) - crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c))) - do s = 1, phase_Nsources(p) - crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c))) - enddo + if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & nonlocalBroken = .true. if(.not. crystallite_todo(g,i,e)) cycle @@ -1567,15 +1543,11 @@ subroutine integrateStateRKCK45 nonlocalBroken = .true. if(.not. crystallite_todo(g,i,e)) exit - call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_todo(g,i,e) = .not. 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) - crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c))) - do s = 1, phase_Nsources(p) - crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c))) - enddo if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & nonlocalBroken = .true. if(.not. crystallite_todo(g,i,e)) exit From 6b11d438424ccc74241719434de678d1d156758d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 08:02:08 +0200 Subject: [PATCH 03/23] handle error checking centrally --- src/constitutive.f90 | 22 +++++++++++++++------- src/crystallite.f90 | 13 ++----------- 2 files changed, 17 insertions(+), 18 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 4402820e8..582222dae 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -327,7 +327,7 @@ module constitutive constitutive_initialFi, & constitutive_SandItsTangents, & constitutive_collectDotState, & - constitutive_collectDeltaState, & + constitutive_deltaState, & constitutive_results contains @@ -794,7 +794,7 @@ end function constitutive_collectDotState !> @brief for constitutive models having an instantaneous change of state !> will return false if delta state is not needed/supported by the constitutive model !-------------------------------------------------------------------------------------------------- -subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el) +function constitutive_deltaState(S, Fe, Fi, ipc, ip, el) result(broken) integer, intent(in) :: & ipc, & !< component-ID of integration point @@ -808,13 +808,18 @@ subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el) Mp integer :: & i, & - instance, of + instance, of, & + phase + logical :: & + broken + Mp = matmul(matmul(transpose(Fi),Fi),S) of = material_phasememberAt(ipc,ip,el) + phase = material_phaseAt(ipc,el) instance = phase_plasticityInstance(material_phaseAt(ipc,el)) - plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) + plasticityType: select case (phase_plasticity(phase)) case (PLASTICITY_KINEHARDENING_ID) plasticityType call plastic_kinehardening_deltaState(Mp,instance,of) @@ -823,10 +828,11 @@ subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el) call plastic_nonlocal_deltaState(Mp,instance,of,ip,el) end select plasticityType + broken = any(IEEE_is_NaN(plasticState(phase)%deltaState(:,of))) - sourceLoop: do i = 1, phase_Nsources(material_phaseAt(ipc,el)) + sourceLoop: do i = 1, phase_Nsources(phase) - sourceType: select case (phase_source(i,material_phaseAt(ipc,el))) + sourceType: select case (phase_source(i,phase)) case (SOURCE_damage_isoBrittle_ID) sourceType call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(ipc,ip,el), Fe, & @@ -834,9 +840,11 @@ subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el) end select sourceType + broken = broken .or. any(IEEE_is_NaN(sourceState(phase)%p(i)%deltaState(:,of))) + enddo SourceLoop -end subroutine constitutive_collectDeltaState +end function constitutive_deltaState !-------------------------------------------------------------------------------------------------- diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 03de4dea3..81119f2e2 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1654,18 +1654,15 @@ logical function stateJump(ipc,ip,el) c = material_phaseMemberAt(ipc,ip,el) p = material_phaseAt(ipc,el) - call constitutive_collectDeltaState(crystallite_S(1:3,1:3,ipc,ip,el), & + stateJump = .not. constitutive_deltaState(crystallite_S(1:3,1:3,ipc,ip,el), & crystallite_Fe(1:3,1:3,ipc,ip,el), & crystallite_Fi(1:3,1:3,ipc,ip,el), & ipc,ip,el) + if(.not. stateJump) return myOffset = plasticState(p)%offsetDeltaState mySize = plasticState(p)%sizeDeltaState - if( any(IEEE_is_NaN(plasticState(p)%deltaState(1:mySize,c)))) then - stateJump = .false. - return - endif plasticState(p)%state(myOffset + 1:myOffset + mySize,c) = & plasticState(p)%state(myOffset + 1:myOffset + mySize,c) + plasticState(p)%deltaState(1:mySize,c) @@ -1673,16 +1670,10 @@ logical function stateJump(ipc,ip,el) do mySource = 1, phase_Nsources(p) myOffset = sourceState(p)%p(mySource)%offsetDeltaState mySize = sourceState(p)%p(mySource)%sizeDeltaState - if (any(IEEE_is_NaN(sourceState(p)%p(mySource)%deltaState(1:mySize,c)))) then - stateJump = .false. - return - endif sourceState(p)%p(mySource)%state(myOffset + 1: myOffset + mySize,c) = & sourceState(p)%p(mySource)%state(myOffset + 1: myOffset + mySize,c) + sourceState(p)%p(mySource)%deltaState(1:mySize,c) enddo - stateJump = .true. - end function stateJump From ec53e4c318328f5d928f45650e09842bda7851ae Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 08:39:49 +0200 Subject: [PATCH 04/23] avoid writing to public variable crystallite_todo --- src/crystallite.f90 | 152 ++++++++++++++++++++++---------------------- 1 file changed, 76 insertions(+), 76 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 81119f2e2..caedf7ef5 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1016,10 +1016,10 @@ subroutine integrateStateFPI real(pReal), dimension(:), allocatable :: plastic_dotState_p1, plastic_dotState_p2 real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState logical :: & - nonlocalBroken + nonlocalBroken, broken nonlocalBroken = .false. - !$OMP PARALLEL DO PRIVATE(sizeDotState,r,zeta,p,c,plastic_dotState_p1, plastic_dotState_p2,source_dotState) + !$OMP PARALLEL DO PRIVATE(sizeDotState,r,zeta,p,c,plastic_dotState_p1, plastic_dotState_p2,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)) @@ -1027,15 +1027,15 @@ subroutine integrateStateFPI p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) - crystallite_todo(g,i,e) = .not. constitutive_collectDotState(crystallite_S(1:3,1:3,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) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) cycle + if(broken) cycle sizeDotState = plasticState(p)%sizeDotState plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & @@ -1060,16 +1060,16 @@ subroutine integrateStateFPI source_dotState(1:sizeDotState,1,s) = sourceState(p)%p(s)%dotState(:,c) enddo - crystallite_todo(g,i,e) = integrateStress(g,i,e) - if(.not. crystallite_todo(g,i,e)) exit iteration + broken = .not. integrateStress(g,i,e) + if(broken) exit iteration - crystallite_todo(g,i,e) = .not. constitutive_collectDotState(crystallite_S(1:3,1:3,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) - if(.not. crystallite_todo(g,i,e)) exit iteration + if(broken) exit iteration sizeDotState = plasticState(p)%sizeDotState zeta = damper(plasticState(p)%dotState(:,c),plastic_dotState_p1,plastic_dotState_p2) @@ -1102,12 +1102,12 @@ subroutine integrateStateFPI enddo if(crystallite_converged(g,i,e)) then - crystallite_todo(g,i,e) = stateJump(g,i,e) + broken = .not. stateJump(g,i,e) exit iteration endif enddo iteration - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. endif @@ -1155,10 +1155,10 @@ subroutine integrateStateEuler s, & sizeDotState logical :: & - nonlocalBroken + nonlocalBroken, broken nonlocalBroken = .false. - !$OMP PARALLEL DO PRIVATE (sizeDotState,p,c) + !$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)) @@ -1166,15 +1166,15 @@ subroutine integrateStateEuler p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) - crystallite_todo(g,i,e) = .not. constitutive_collectDotState(crystallite_S(1:3,1:3,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) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) cycle + if(broken) cycle sizeDotState = plasticState(p)%sizeDotState plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & @@ -1187,16 +1187,16 @@ subroutine integrateStateEuler * crystallite_subdt(g,i,e) enddo - crystallite_todo(g,i,e) = stateJump(g,i,e) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + broken = .not. stateJump(g,i,e) + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) cycle + if(broken) cycle - crystallite_todo(g,i,e) = integrateStress(g,i,e) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + broken = .not. integrateStress(g,i,e) + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - crystallite_converged(g,i,e) = crystallite_todo(g,i,e) + crystallite_converged(g,i,e) = .not. broken endif enddo; enddo; enddo @@ -1221,13 +1221,13 @@ subroutine integrateStateAdaptiveEuler s, & sizeDotState logical :: & - nonlocalBroken + nonlocalBroken, 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) + !$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)) @@ -1235,15 +1235,15 @@ subroutine integrateStateAdaptiveEuler p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) - crystallite_todo(g,i,e) = .not. constitutive_collectDotState(crystallite_S(1:3,1:3,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) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) cycle + if(broken) cycle sizeDotState = plasticState(p)%sizeDotState @@ -1259,25 +1259,25 @@ subroutine integrateStateAdaptiveEuler + sourceState(p)%p(s)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) enddo - crystallite_todo(g,i,e) = stateJump(g,i,e) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + broken = .not. stateJump(g,i,e) + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) cycle + if(broken) cycle - crystallite_todo(g,i,e) = integrateStress(g,i,e) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + broken = .not. integrateStress(g,i,e) + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) cycle + if(broken) cycle - crystallite_todo(g,i,e) = .not. constitutive_collectDotState(crystallite_S(1:3,1:3,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) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) cycle + if(broken) cycle sizeDotState = plasticState(p)%sizeDotState @@ -1333,12 +1333,12 @@ subroutine integrateStateRK4 s, & sizeDotState logical :: & - nonlocalBroken + nonlocalBroken, broken real(pReal), dimension(constitutive_source_maxSizeDotState,4,maxval(phase_Nsources)) :: source_RK4dotState real(pReal), dimension(constitutive_plasticity_maxSizeDotState,4) :: plastic_RK4dotState nonlocalBroken = .false. - !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,source_RK4dotState,plastic_RK4dotState) + !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,source_RK4dotState,plastic_RK4dotState,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)) @@ -1346,15 +1346,15 @@ subroutine integrateStateRK4 p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) - crystallite_todo(g,i,e) = .not. constitutive_collectDotState(crystallite_S(1:3,1:3,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) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) cycle + if(broken) cycle do stage = 1,3 sizeDotState = plasticState(p)%sizeDotState @@ -1388,24 +1388,24 @@ subroutine integrateStateRK4 * crystallite_subdt(g,i,e) enddo - crystallite_todo(g,i,e) = integrateStress(g,i,e,CC(stage)) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + broken = .not. integrateStress(g,i,e,CC(stage)) + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) exit + if(broken) exit - crystallite_todo(g,i,e) = .not. constitutive_collectDotState(crystallite_S(1:3,1:3,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)*CC(stage), g,i,e) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) exit + if(broken) exit enddo - if(.not. crystallite_todo(g,i,e)) cycle + if(broken) cycle sizeDotState = plasticState(p)%sizeDotState @@ -1427,15 +1427,15 @@ subroutine integrateStateRK4 * crystallite_subdt(g,i,e) enddo - crystallite_todo(g,i,e) = stateJump(g,i,e) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + broken = .not. stateJump(g,i,e) + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) cycle + if(broken) cycle - crystallite_todo(g,i,e) = integrateStress(g,i,e) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + broken = .not. integrateStress(g,i,e) + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - crystallite_converged(g,i,e) = crystallite_todo(g,i,e) ! consider converged if not broken + crystallite_converged(g,i,e) = .not. broken endif enddo; enddo; enddo @@ -1483,12 +1483,12 @@ subroutine integrateStateRKCK45 s, & sizeDotState logical :: & - nonlocalBroken + nonlocalBroken, broken real(pReal), dimension(constitutive_source_maxSizeDotState,6,maxval(phase_Nsources)) :: source_RKdotState real(pReal), dimension(constitutive_plasticity_maxSizeDotState,6) :: plastic_RKdotState nonlocalBroken = .false. - !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,plastic_RKdotState,source_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)) @@ -1496,15 +1496,15 @@ subroutine integrateStateRKCK45 p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) - crystallite_todo(g,i,e) = .not. constitutive_collectDotState(crystallite_S(1:3,1:3,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) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) cycle + if(broken) cycle do stage = 1,5 sizeDotState = plasticState(p)%sizeDotState @@ -1538,23 +1538,23 @@ subroutine integrateStateRKCK45 * crystallite_subdt(g,i,e) enddo - crystallite_todo(g,i,e) = integrateStress(g,i,e,CC(stage)) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + broken = .not. integrateStress(g,i,e,CC(stage)) + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) exit + if(broken) exit - crystallite_todo(g,i,e) = .not. constitutive_collectDotState(crystallite_S(1:3,1:3,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)*CC(stage), g,i,e) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) exit + if(broken) exit enddo - if(.not. crystallite_todo(g,i,e)) cycle + if(broken) cycle sizeDotState = plasticState(p)%sizeDotState @@ -1563,7 +1563,7 @@ subroutine integrateStateRKCK45 plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & + plasticState(p)%dotState (1:sizeDotState,c) & * crystallite_subdt(g,i,e) - crystallite_todo(g,i,e) = converged( matmul(plastic_RKdotState(1:sizeDotState,1:6),DB) & + broken = .not. converged( matmul(plastic_RKdotState(1:sizeDotState,1:6),DB) & * crystallite_subdt(g,i,e), & plasticState(p)%state(1:sizeDotState,c), & plasticState(p)%atol(1:sizeDotState)) @@ -1576,25 +1576,25 @@ subroutine integrateStateRKCK45 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) - crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. & + broken = broken .and. .not. & converged(matmul(source_RKdotState(1:sizeDotState,1:6,s),DB) & * crystallite_subdt(g,i,e), & sourceState(p)%p(s)%state(1:sizeDotState,c), & sourceState(p)%p(s)%atol(1:sizeDotState)) enddo - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) cycle + if(broken) cycle - crystallite_todo(g,i,e) = stateJump(g,i,e) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + broken = .not. stateJump(g,i,e) + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) cycle + if(broken) cycle - crystallite_todo(g,i,e) = integrateStress(g,i,e) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + broken = .not. integrateStress(g,i,e) + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - crystallite_converged(g,i,e) = crystallite_todo(g,i,e) ! consider converged if not broken + crystallite_converged(g,i,e) = .not. broken endif enddo; enddo; enddo From d50d55cef3ac57d0e86386c8103be832f889b8c2 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 08:54:20 +0200 Subject: [PATCH 05/23] avoid public variables --- src/crystallite.f90 | 50 ++++++++++++++++++++++++++------------------- 1 file changed, 29 insertions(+), 21 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index caedf7ef5..908d57ef8 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -71,7 +71,6 @@ module crystallite crystallite_requested !< used by upper level (homogenization) to request crystallite calculation logical, dimension(:,:,:), allocatable :: & crystallite_converged, & !< convergence flag - crystallite_todo, & !< flag to indicate need for further computation crystallite_localPlasticity !< indicates this grain to have purely local constitutive law type :: tOutput !< new requested output (per phase) @@ -98,7 +97,7 @@ module crystallite type(tNumerics) :: num ! numerics parameters. Better name? - procedure(), pointer :: integrateState + procedure(integrateStateFPI), pointer :: integrateState public :: & crystallite_init, & @@ -161,7 +160,6 @@ subroutine crystallite_init allocate(crystallite_localPlasticity(cMax,iMax,eMax), source=.true.) allocate(crystallite_requested(cMax,iMax,eMax), source=.false.) - allocate(crystallite_todo(cMax,iMax,eMax), source=.false.) allocate(crystallite_converged(cMax,iMax,eMax), source=.true.) num%subStepMinCryst = config_numerics%getFloat('substepmincryst', defaultVal=1.0e-3_pReal) @@ -301,6 +299,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) e, & !< counter in element loop startIP, endIP, & s + logical, dimension(homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: todo !ToDo: need to set some values to false in hase of different Ngrains #ifdef DEBUG if (iand(debug_level(debug_crystallite),debug_levelSelective) /= 0 & @@ -344,7 +343,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) 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 - crystallite_todo(c,i,e) = .true. + todo(c,i,e) = .true. crystallite_converged(c,i,e) = .false. ! pretend failed step of 1/subStepSizeCryst endif homogenizationRequestsCalculation enddo; enddo @@ -361,7 +360,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) endif singleRun NiterationCrystallite = 0 - cutbackLooping: do while (any(crystallite_todo(:,startIP:endIP,FEsolving_execELem(1):FEsolving_execElem(2)))) + cutbackLooping: do while (any(todo(:,startIP:endIP,FEsolving_execELem(1):FEsolving_execElem(2)))) NiterationCrystallite = NiterationCrystallite + 1 #ifdef DEBUG @@ -380,8 +379,8 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) crystallite_subStep(c,i,e) = min(1.0_pReal - crystallite_subFrac(c,i,e), & num%stepIncreaseCryst * crystallite_subStep(c,i,e)) - crystallite_todo(c,i,e) = crystallite_subStep(c,i,e) > 0.0_pReal ! still time left to integrate on? - if (crystallite_todo(c,i,e)) then + 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) @@ -415,12 +414,12 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) enddo ! cant restore dotState here, since not yet calculated in first cutback after initialization - crystallite_todo(c,i,e) = crystallite_subStep(c,i,e) > num%subStepMinCryst ! still on track or already done (beyond repair) + todo(c,i,e) = crystallite_subStep(c,i,e) > num%subStepMinCryst ! still on track or already done (beyond repair) endif !-------------------------------------------------------------------------------------------------- ! prepare for integration - if (crystallite_todo(c,i,e)) then + if (todo(c,i,e)) then crystallite_subF(1:3,1:3,c,i,e) = crystallite_subF0(1:3,1:3,c,i,e) & + crystallite_subStep(c,i,e) *( crystallite_partionedF (1:3,1:3,c,i,e) & -crystallite_partionedF0(1:3,1:3,c,i,e)) @@ -438,9 +437,9 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) !-------------------------------------------------------------------------------------------------- ! integrate --- requires fully defined state array (basic + dependent state) - if (any(crystallite_todo)) call integrateState ! TODO: unroll into proper elementloop to avoid N^2 for single point evaluation + 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 - crystallite_todo = .true. ! TODO: again unroll this into proper elementloop to avoid N^2 for single point evaluation + todo = .true. ! TODO: again unroll this into proper elementloop to avoid N^2 for single point evaluation enddo cutbackLooping @@ -998,8 +997,9 @@ end function integrateStress !> @brief integrate stress, state with adaptive 1st order explicit Euler method !> using Fixed Point Iteration to adapt the stepsize !-------------------------------------------------------------------------------------------------- -subroutine integrateStateFPI +subroutine integrateStateFPI(todo) + logical, dimension(:,:,:), intent(in) :: todo integer :: & NiterationState, & !< number of iterations in state loop e, & !< element index in element loop @@ -1023,7 +1023,7 @@ subroutine integrateStateFPI 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(crystallite_todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then + if(todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) @@ -1144,7 +1144,9 @@ end subroutine integrateStateFPI !-------------------------------------------------------------------------------------------------- !> @brief integrate state with 1st order explicit Euler method !-------------------------------------------------------------------------------------------------- -subroutine integrateStateEuler +subroutine integrateStateEuler(todo) + + logical, dimension(:,:,:), intent(in) :: todo integer :: & e, & !< element index in element loop @@ -1162,7 +1164,7 @@ subroutine integrateStateEuler 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(crystallite_todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then + if(todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) @@ -1210,7 +1212,9 @@ end subroutine integrateStateEuler !-------------------------------------------------------------------------------------------------- !> @brief integrate stress, state with 1st order Euler method with adaptive step size !-------------------------------------------------------------------------------------------------- -subroutine integrateStateAdaptiveEuler +subroutine integrateStateAdaptiveEuler(todo) + + logical, dimension(:,:,:), intent(in) :: todo integer :: & e, & ! element index in element loop @@ -1231,7 +1235,7 @@ subroutine integrateStateAdaptiveEuler 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(crystallite_todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then + if(todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) @@ -1309,7 +1313,9 @@ end subroutine integrateStateAdaptiveEuler !-------------------------------------------------------------------------------------------------- !> @brief integrate stress, state with 4th order explicit Runge Kutta method !-------------------------------------------------------------------------------------------------- -subroutine integrateStateRK4 +subroutine integrateStateRK4(todo) + + logical, dimension(:,:,:), intent(in) :: todo real(pReal), dimension(3,3), parameter :: & A = reshape([& @@ -1342,7 +1348,7 @@ subroutine integrateStateRK4 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(crystallite_todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then + if(todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) @@ -1450,7 +1456,9 @@ end subroutine integrateStateRK4 !> @brief integrate stress, state with 5th order Runge-Kutta Cash-Karp method with !> adaptive step size (use 5th order solution to advance = "local extrapolation") !-------------------------------------------------------------------------------------------------- -subroutine integrateStateRKCK45 +subroutine integrateStateRKCK45(todo) + + logical, dimension(:,:,:), intent(in) :: todo real(pReal), dimension(5,5), parameter :: & A = reshape([& @@ -1492,7 +1500,7 @@ subroutine integrateStateRKCK45 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(crystallite_todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then + if(todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) From 6eee8f34ac830f4d2db934eb92db66211a4eb886 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 09:01:03 +0200 Subject: [PATCH 06/23] homogeneous mesh --- src/crystallite.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 908d57ef8..4e475343b 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -269,7 +269,7 @@ subroutine crystallite_init #ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0) then write(6,'(a42,1x,i10)') ' # of elements: ', eMax - write(6,'(a42,1x,i10)') 'max # of integration points/element: ', iMax + write(6,'(a42,1x,i10)') ' # of integration points/element: ', iMax write(6,'(a42,1x,i10)') 'max # of constituents/integration point: ', cMax write(6,'(a42,1x,i10)') ' # of nonlocal constituents: ',count(.not. crystallite_localPlasticity) flush(6) From 079596f7cdcaa9ac69c7ccaa5260c51bb1725c30 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 09:23:43 +0200 Subject: [PATCH 07/23] unix standard 0/.false. = OK, 1/.true. = not OK --- src/crystallite.f90 | 46 ++++++++++++++++++++++----------------------- 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 4e475343b..4fe8c75dd 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -776,7 +776,7 @@ end subroutine crystallite_results !> @brief calculation of stress (P) with time integration based on a residuum in Lp and !> intermediate acceleration of the Newton-Raphson correction !-------------------------------------------------------------------------------------------------- -logical function integrateStress(ipc,ip,el,timeFraction) +function integrateStress(ipc,ip,el,timeFraction) result(broken) integer, intent(in):: el, & ! element index ip, & ! integration point index @@ -833,11 +833,11 @@ logical function integrateStress(ipc,ip,el,timeFraction) p, & jacoCounterLp, & jacoCounterLi ! counters to check for Jacobian update - logical :: error + logical :: error,broken external :: & dgesv - integrateStress = .false. + broken = .true. if (present(timeFraction)) then dt = crystallite_subdt(ipc,ip,el) * timeFraction @@ -981,7 +981,6 @@ logical function integrateStress(ipc,ip,el,timeFraction) call math_invert33(Fp_new,devNull,error,invFp_new) if (error) return ! error - integrateStress = .true. crystallite_P (1:3,1:3,ipc,ip,el) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new))) crystallite_S (1:3,1:3,ipc,ip,el) = S crystallite_Lp (1:3,1:3,ipc,ip,el) = Lpguess @@ -989,6 +988,7 @@ logical function integrateStress(ipc,ip,el,timeFraction) crystallite_Fp (1:3,1:3,ipc,ip,el) = Fp_new / math_det33(Fp_new)**(1.0_pReal/3.0_pReal) ! regularize crystallite_Fi (1:3,1:3,ipc,ip,el) = Fi_new crystallite_Fe (1:3,1:3,ipc,ip,el) = matmul(matmul(F,invFp_new),invFi_new) + broken = .false. end function integrateStress @@ -1060,7 +1060,7 @@ subroutine integrateStateFPI(todo) source_dotState(1:sizeDotState,1,s) = sourceState(p)%p(s)%dotState(:,c) enddo - broken = .not. integrateStress(g,i,e) + broken = integrateStress(g,i,e) if(broken) exit iteration broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & @@ -1068,7 +1068,6 @@ subroutine integrateStateFPI(todo) crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e) - if(broken) exit iteration sizeDotState = plasticState(p)%sizeDotState @@ -1102,7 +1101,7 @@ subroutine integrateStateFPI(todo) enddo if(crystallite_converged(g,i,e)) then - broken = .not. stateJump(g,i,e) + broken = stateJump(g,i,e) exit iteration endif @@ -1189,12 +1188,12 @@ subroutine integrateStateEuler(todo) * crystallite_subdt(g,i,e) enddo - broken = .not. stateJump(g,i,e) + broken = stateJump(g,i,e) if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. if(broken) cycle - broken = .not. integrateStress(g,i,e) + broken = integrateStress(g,i,e) if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. @@ -1263,12 +1262,12 @@ subroutine integrateStateAdaptiveEuler(todo) + sourceState(p)%p(s)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) enddo - broken = .not. stateJump(g,i,e) + broken = stateJump(g,i,e) if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. if(broken) cycle - broken = .not. integrateStress(g,i,e) + broken = integrateStress(g,i,e) if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. if(broken) cycle @@ -1394,7 +1393,7 @@ subroutine integrateStateRK4(todo) * crystallite_subdt(g,i,e) enddo - broken = .not. integrateStress(g,i,e,CC(stage)) + broken = integrateStress(g,i,e,CC(stage)) if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. if(broken) exit @@ -1433,12 +1432,12 @@ subroutine integrateStateRK4(todo) * crystallite_subdt(g,i,e) enddo - broken = .not. stateJump(g,i,e) + broken = stateJump(g,i,e) if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. if(broken) cycle - broken = .not. integrateStress(g,i,e) + broken = integrateStress(g,i,e) if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. crystallite_converged(g,i,e) = .not. broken @@ -1546,7 +1545,7 @@ subroutine integrateStateRKCK45(todo) * crystallite_subdt(g,i,e) enddo - broken = .not. integrateStress(g,i,e,CC(stage)) + broken = integrateStress(g,i,e,CC(stage)) if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. if(broken) exit @@ -1594,12 +1593,12 @@ subroutine integrateStateRKCK45(todo) nonlocalBroken = .true. if(broken) cycle - broken = .not. stateJump(g,i,e) + broken = stateJump(g,i,e) if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. if(broken) cycle - broken = .not. integrateStress(g,i,e) + broken = integrateStress(g,i,e) if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. crystallite_converged(g,i,e) = .not. broken @@ -1645,7 +1644,7 @@ end function converged !> @brief calculates a jump in the state according to the current state and the current stress !> returns true, if state jump was successfull or not needed. false indicates NaN in delta state !-------------------------------------------------------------------------------------------------- -logical function stateJump(ipc,ip,el) +function stateJump(ipc,ip,el) result(broken) integer, intent(in):: & el, & ! element index @@ -1658,15 +1657,16 @@ logical function stateJump(ipc,ip,el) mySource, & myOffset, & mySize + logical :: broken c = material_phaseMemberAt(ipc,ip,el) p = material_phaseAt(ipc,el) - stateJump = .not. constitutive_deltaState(crystallite_S(1:3,1:3,ipc,ip,el), & - crystallite_Fe(1:3,1:3,ipc,ip,el), & - crystallite_Fi(1:3,1:3,ipc,ip,el), & - ipc,ip,el) - if(.not. stateJump) return + broken = constitutive_deltaState(crystallite_S(1:3,1:3,ipc,ip,el), & + crystallite_Fe(1:3,1:3,ipc,ip,el), & + crystallite_Fi(1:3,1:3,ipc,ip,el), & + ipc,ip,el) + if(broken) return myOffset = plasticState(p)%offsetDeltaState mySize = plasticState(p)%sizeDeltaState From ce61606c0bd2563b0adeaa10d6d03017b290a5f7 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 09:34:49 +0200 Subject: [PATCH 08/23] not needed --- src/material.f90 | 1 - 1 file changed, 1 deletion(-) diff --git a/src/material.f90 b/src/material.f90 index aefe44878..2a447c036 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -11,7 +11,6 @@ module material use results use IO use debug - use numerics use rotations use discretization From b996b6c42ec1e0954b4e4896cd7950b99f0890e5 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 09:56:59 +0200 Subject: [PATCH 09/23] merge functionality --- src/constitutive_plastic_disloUCLA.f90 | 2 +- src/constitutive_plastic_dislotwin.f90 | 2 +- src/constitutive_plastic_isotropic.f90 | 2 +- src/constitutive_plastic_kinehardening.f90 | 2 +- src/constitutive_plastic_none.f90 | 2 +- src/constitutive_plastic_nonlocal.f90 | 2 +- src/constitutive_plastic_phenopowerlaw.f90 | 2 +- src/material.f90 | 65 ++++++---------------- src/source_damage_anisoBrittle.f90 | 2 +- src/source_damage_anisoDuctile.f90 | 2 +- src/source_damage_isoBrittle.f90 | 2 +- src/source_damage_isoDuctile.f90 | 2 +- src/source_thermal_dissipation.f90 | 2 +- src/source_thermal_externalheat.f90 | 2 +- 14 files changed, 31 insertions(+), 60 deletions(-) diff --git a/src/constitutive_plastic_disloUCLA.f90 b/src/constitutive_plastic_disloUCLA.f90 index 6be86f266..90a933910 100644 --- a/src/constitutive_plastic_disloUCLA.f90 +++ b/src/constitutive_plastic_disloUCLA.f90 @@ -209,7 +209,7 @@ module subroutine plastic_disloUCLA_init sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl sizeState = sizeDotState - call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0) + call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,0) !-------------------------------------------------------------------------------------------------- ! state aliases and initialization diff --git a/src/constitutive_plastic_dislotwin.f90 b/src/constitutive_plastic_dislotwin.f90 index d36e08846..abc46a45e 100644 --- a/src/constitutive_plastic_dislotwin.f90 +++ b/src/constitutive_plastic_dislotwin.f90 @@ -399,7 +399,7 @@ module subroutine plastic_dislotwin_init + size(['f_tr']) * prm%sum_N_tr sizeState = sizeDotState - call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0) + call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,0) !-------------------------------------------------------------------------------------------------- ! locally defined state aliases and initialization of state0 and atol diff --git a/src/constitutive_plastic_isotropic.f90 b/src/constitutive_plastic_isotropic.f90 index 94fc9817d..ecf029124 100644 --- a/src/constitutive_plastic_isotropic.f90 +++ b/src/constitutive_plastic_isotropic.f90 @@ -117,7 +117,7 @@ module subroutine plastic_isotropic_init sizeDotState = size(['xi ','accumulated_shear']) sizeState = sizeDotState - call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0) + call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,0) !-------------------------------------------------------------------------------------------------- ! state aliases and initialization diff --git a/src/constitutive_plastic_kinehardening.f90 b/src/constitutive_plastic_kinehardening.f90 index 5843f5b5e..36b1eedf9 100644 --- a/src/constitutive_plastic_kinehardening.f90 +++ b/src/constitutive_plastic_kinehardening.f90 @@ -164,7 +164,7 @@ module subroutine plastic_kinehardening_init sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl sizeState = sizeDotState + sizeDeltaState - call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,sizeDeltaState) + call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,sizeDeltaState) !-------------------------------------------------------------------------------------------------- ! state aliases and initialization diff --git a/src/constitutive_plastic_none.f90 b/src/constitutive_plastic_none.f90 index 7ff1c76f7..667fe5638 100644 --- a/src/constitutive_plastic_none.f90 +++ b/src/constitutive_plastic_none.f90 @@ -29,7 +29,7 @@ module subroutine plastic_none_init if (phase_plasticity(p) /= PLASTICITY_NONE_ID) cycle NipcMyPhase = count(material_phaseAt == p) * discretization_nIP - call material_allocatePlasticState(p,NipcMyPhase,0,0,0) + call material_allocateState(plasticState(p),NipcMyPhase,0,0,0) enddo diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index 58218ac00..35cc7bdd1 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -384,7 +384,7 @@ module subroutine plastic_nonlocal_init 'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%sum_N_sl !< other dependent state variables that are not updated by microstructure sizeDeltaState = sizeDotState - call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,sizeDeltaState) + call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,sizeDeltaState) plasticState(p)%nonlocal = .true. plasticState(p)%offsetDeltaState = 0 ! ToDo: state structure does not follow convention diff --git a/src/constitutive_plastic_phenopowerlaw.f90 b/src/constitutive_plastic_phenopowerlaw.f90 index a980d6106..12a30478a 100644 --- a/src/constitutive_plastic_phenopowerlaw.f90 +++ b/src/constitutive_plastic_phenopowerlaw.f90 @@ -213,7 +213,7 @@ module subroutine plastic_phenopowerlaw_init + size(['xi_tw ','gamma_tw']) * prm%sum_N_tw sizeState = sizeDotState - call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0) + call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,0) !-------------------------------------------------------------------------------------------------- ! state aliases and initialization diff --git a/src/material.f90 b/src/material.f90 index 2a447c036..749c9a3d8 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -173,8 +173,7 @@ module material public :: & material_init, & - material_allocatePlasticState, & - material_allocateSourceState, & + material_allocateState, & ELASTICITY_HOOKE_ID ,& PLASTICITY_NONE_ID, & PLASTICITY_ISOTROPIC_ID, & @@ -699,63 +698,35 @@ end subroutine material_parseTexture !-------------------------------------------------------------------------------------------------- -!> @brief allocates the plastic state of a phase +!> @brief Allocate the components of the state structure for a given phase !-------------------------------------------------------------------------------------------------- -subroutine material_allocatePlasticState(phase,NipcMyPhase,& - sizeState,sizeDotState,sizeDeltaState) +subroutine material_allocateState(state, & + NipcMyPhase,sizeState,sizeDotState,sizeDeltaState) + class(tState), intent(out) :: & + state integer, intent(in) :: & - phase, & NipcMyPhase, & sizeState, & sizeDotState, & sizeDeltaState - plasticState(phase)%sizeState = sizeState - plasticState(phase)%sizeDotState = sizeDotState - plasticState(phase)%sizeDeltaState = sizeDeltaState - plasticState(phase)%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition + state%sizeState = sizeState + state%sizeDotState = sizeDotState + state%sizeDeltaState = sizeDeltaState + state%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition - allocate(plasticState(phase)%atol (sizeState), source=0.0_pReal) - allocate(plasticState(phase)%state0 (sizeState,NipcMyPhase), source=0.0_pReal) - allocate(plasticState(phase)%partionedState0 (sizeState,NipcMyPhase), source=0.0_pReal) - allocate(plasticState(phase)%subState0 (sizeState,NipcMyPhase), source=0.0_pReal) - allocate(plasticState(phase)%state (sizeState,NipcMyPhase), source=0.0_pReal) + allocate(state%atol (sizeState), source=0.0_pReal) + allocate(state%state0 (sizeState,NipcMyPhase), source=0.0_pReal) + allocate(state%partionedState0(sizeState,NipcMyPhase), source=0.0_pReal) + allocate(state%subState0 (sizeState,NipcMyPhase), source=0.0_pReal) + allocate(state%state (sizeState,NipcMyPhase), source=0.0_pReal) - allocate(plasticState(phase)%dotState (sizeDotState,NipcMyPhase),source=0.0_pReal) + allocate(state%dotState (sizeDotState,NipcMyPhase), source=0.0_pReal) - allocate(plasticState(phase)%deltaState (sizeDeltaState,NipcMyPhase),source=0.0_pReal) + allocate(state%deltaState(sizeDeltaState,NipcMyPhase), source=0.0_pReal) -end subroutine material_allocatePlasticState +end subroutine material_allocateState -!-------------------------------------------------------------------------------------------------- -!> @brief allocates the source state of a phase -!-------------------------------------------------------------------------------------------------- -subroutine material_allocateSourceState(phase,of,NipcMyPhase,& - sizeState,sizeDotState,sizeDeltaState) - - integer, intent(in) :: & - phase, & - of, & - NipcMyPhase, & - sizeState, sizeDotState,sizeDeltaState - - sourceState(phase)%p(of)%sizeState = sizeState - sourceState(phase)%p(of)%sizeDotState = sizeDotState - sourceState(phase)%p(of)%sizeDeltaState = sizeDeltaState - sourceState(phase)%p(of)%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition - - allocate(sourceState(phase)%p(of)%atol (sizeState), source=0.0_pReal) - allocate(sourceState(phase)%p(of)%state0 (sizeState,NipcMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(of)%partionedState0 (sizeState,NipcMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(of)%subState0 (sizeState,NipcMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(of)%state (sizeState,NipcMyPhase), source=0.0_pReal) - - allocate(sourceState(phase)%p(of)%dotState (sizeDotState,NipcMyPhase),source=0.0_pReal) - - allocate(sourceState(phase)%p(of)%deltaState (sizeDeltaState,NipcMyPhase),source=0.0_pReal) - -end subroutine material_allocateSourceState - end module material diff --git a/src/source_damage_anisoBrittle.f90 b/src/source_damage_anisoBrittle.f90 index 3978be959..b3af24f38 100644 --- a/src/source_damage_anisoBrittle.f90 +++ b/src/source_damage_anisoBrittle.f90 @@ -107,7 +107,7 @@ subroutine source_damage_anisoBrittle_init if (any(prm%critDisp < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_critDisp' NipcMyPhase = count(material_phaseAt==p) * discretization_nIP - call material_allocateSourceState(p,sourceOffset,NipcMyPhase,1,1,0) + call material_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0) sourceState(p)%p(sourceOffset)%atol = config%getFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal) if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol' diff --git a/src/source_damage_anisoDuctile.f90 b/src/source_damage_anisoDuctile.f90 index 2b818e2cf..79cc0c2f7 100644 --- a/src/source_damage_anisoDuctile.f90 +++ b/src/source_damage_anisoDuctile.f90 @@ -89,7 +89,7 @@ subroutine source_damage_anisoDuctile_init if (any(prm%critPlasticStrain < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_criticalplasticstrain' NipcMyPhase=count(material_phaseAt==p) * discretization_nIP - call material_allocateSourceState(p,sourceOffset,NipcMyPhase,1,1,0) + call material_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0) sourceState(p)%p(sourceOffset)%atol = config%getFloat('anisoductile_atol',defaultVal=1.0e-3_pReal) if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol' diff --git a/src/source_damage_isoBrittle.f90 b/src/source_damage_isoBrittle.f90 index ed6d89a89..9eacb4516 100644 --- a/src/source_damage_isoBrittle.f90 +++ b/src/source_damage_isoBrittle.f90 @@ -83,7 +83,7 @@ subroutine source_damage_isoBrittle_init if (prm%critStrainEnergy <= 0.0_pReal) extmsg = trim(extmsg)//' isobrittle_criticalstrainenergy' NipcMyPhase = count(material_phaseAt==p) * discretization_nIP - call material_allocateSourceState(p,sourceOffset,NipcMyPhase,1,1,1) + call material_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,1) sourceState(p)%p(sourceOffset)%atol = config%getFloat('isobrittle_atol',defaultVal=1.0e-3_pReal) if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol' diff --git a/src/source_damage_isoDuctile.f90 b/src/source_damage_isoDuctile.f90 index 7024e595a..96754725d 100644 --- a/src/source_damage_isoDuctile.f90 +++ b/src/source_damage_isoDuctile.f90 @@ -82,7 +82,7 @@ subroutine source_damage_isoDuctile_init if (prm%critPlasticStrain <= 0.0_pReal) extmsg = trim(extmsg)//' isoductile_criticalplasticstrain' NipcMyPhase=count(material_phaseAt==p) * discretization_nIP - call material_allocateSourceState(p,sourceOffset,NipcMyPhase,1,1,0) + call material_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0) sourceState(p)%p(sourceOffset)%atol = config%getFloat('isoductile_atol',defaultVal=1.0e-3_pReal) if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol' diff --git a/src/source_thermal_dissipation.f90 b/src/source_thermal_dissipation.f90 index 521c79077..c323e68b5 100644 --- a/src/source_thermal_dissipation.f90 +++ b/src/source_thermal_dissipation.f90 @@ -67,7 +67,7 @@ subroutine source_thermal_dissipation_init prm%kappa = config%getFloat('dissipation_coldworkcoeff') NipcMyPhase = count(material_phaseAt==p) * discretization_nIP - call material_allocateSourceState(p,sourceOffset,NipcMyPhase,0,0,0) + call material_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,0,0,0) end associate enddo diff --git a/src/source_thermal_externalheat.f90 b/src/source_thermal_externalheat.f90 index ade13bef2..06b8a5197 100644 --- a/src/source_thermal_externalheat.f90 +++ b/src/source_thermal_externalheat.f90 @@ -74,7 +74,7 @@ subroutine source_thermal_externalheat_init prm%heat_rate = config%getFloats('externalheat_rate',requiredSize = size(prm%time)) NipcMyPhase = count(material_phaseAt==p) * discretization_nIP - call material_allocateSourceState(p,sourceOffset,NipcMyPhase,1,1,0) + call material_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0) end associate enddo From 79012c9ffb0c586ac916b42b73452acdd1d41849 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 10:00:12 +0200 Subject: [PATCH 10/23] not needed --- src/prec.f90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/prec.f90 b/src/prec.f90 index dd39ddc05..646f7dd69 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -53,8 +53,7 @@ module prec logical :: & nonlocal = .false. real(pReal), pointer, dimension(:,:) :: & - slipRate, & !< slip rate - accumulatedSlip !< accumulated plastic slip + slipRate !< slip rate end type type :: tSourceState From cde8c65bd18b9f41bd669f8db673f86066bcdc74 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 10:57:53 +0200 Subject: [PATCH 11/23] better store data locally --- src/crystallite.f90 | 12 +++++++++--- src/numerics.f90 | 23 ++++++++++------------- 2 files changed, 19 insertions(+), 16 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 4fe8c75dd..7b3ee4115 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -15,7 +15,6 @@ module crystallite use DAMASK_interface use config use debug - use numerics use rotations use math use FEsolving @@ -83,7 +82,8 @@ module crystallite integer :: & iJacoLpresiduum, & !< frequency of Jacobian update of residuum in Lp nState, & !< state loop limit - nStress !< stress loop limit + nStress, & !< stress loop limit + integrator !< integration scheme (ToDo: better use a string) real(pReal) :: & subStepMinCryst, & !< minimum (relative) size of sub-step allowed during cutback subStepSizeCryst, & !< size of first substep when cutback @@ -175,6 +175,8 @@ subroutine crystallite_init num%iJacoLpresiduum = config_numerics%getInt ('ijacolpresiduum', defaultVal=1) + num%integrator = config_numerics%getInt ('integrator', defaultVal=1) + num%nState = config_numerics%getInt ('nstate', defaultVal=20) num%nStress = config_numerics%getInt ('nstress', defaultVal=40) @@ -191,10 +193,14 @@ subroutine crystallite_init if(num%iJacoLpresiduum < 1) call IO_error(301,ext_msg='iJacoLpresiduum') + if(num%integrator < 1 .or. num%integrator > 5) & + call IO_error(301,ext_msg='integrator') + if(num%nState < 1) call IO_error(301,ext_msg='nState') if(num%nStress< 1) call IO_error(301,ext_msg='nStress') - select case(numerics_integrator) + + select case(num%integrator) case(1) integrateState => integrateStateFPI case(2) diff --git a/src/numerics.f90 b/src/numerics.f90 index a29601322..8d242c71d 100644 --- a/src/numerics.f90 +++ b/src/numerics.f90 @@ -20,8 +20,7 @@ module numerics iJacoStiffness = 1, & !< frequency of stiffness update randomSeed = 0, & !< fixed seeding for pseudo-random number generator, Default 0: use random seed worldrank = 0, & !< MPI worldrank (/=0 for MPI simulations only) - worldsize = 1, & !< MPI worldsize (/=1 for MPI simulations only) - numerics_integrator = 1 !< method used for state integration Default 1: fix-point iteration + worldsize = 1 !< MPI worldsize (/=1 for MPI simulations only) integer(4), protected, public :: & DAMASK_NumThreadsInt = 0 !< value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive real(pReal), protected, public :: & @@ -134,8 +133,6 @@ subroutine numerics_init defgradTolerance = IO_floatValue(line,chunkPos,2) case ('ijacostiffness') iJacoStiffness = IO_intValue(line,chunkPos,2) - case ('integrator') - numerics_integrator = IO_intValue(line,chunkPos,2) case ('usepingpong') usepingpong = IO_intValue(line,chunkPos,2) > 0 case ('unitlength') @@ -176,6 +173,11 @@ subroutine numerics_init case ('maxstaggerediter') stagItMax = IO_intValue(line,chunkPos,2) +#ifdef PETSC + case ('petsc_options') + petsc_options = trim(line(chunkPos(4):)) +#endif + !-------------------------------------------------------------------------------------------------- ! spectral parameters #ifdef Grid @@ -187,8 +189,6 @@ subroutine numerics_init err_stress_tolrel = IO_floatValue(line,chunkPos,2) case ('err_stress_tolabs') err_stress_tolabs = IO_floatValue(line,chunkPos,2) - case ('petsc_options') - petsc_options = trim(line(chunkPos(4):)) case ('err_curl_tolabs') err_curl_tolAbs = IO_floatValue(line,chunkPos,2) case ('err_curl_tolrel') @@ -206,8 +206,6 @@ subroutine numerics_init integrationorder = IO_intValue(line,chunkPos,2) case ('structorder') structorder = IO_intValue(line,chunkPos,2) - case ('petsc_options') - petsc_options = trim(line(chunkPos(4):)) case ('bbarstabilisation') BBarStabilisation = IO_intValue(line,chunkPos,2) > 0 #endif @@ -223,7 +221,6 @@ subroutine numerics_init ! writing parameters to output write(6,'(a24,1x,es8.1)') ' defgradTolerance: ',defgradTolerance write(6,'(a24,1x,i8)') ' iJacoStiffness: ',iJacoStiffness - write(6,'(a24,1x,i8)') ' integrator: ',numerics_integrator write(6,'(a24,1x,L8)') ' use ping pong scheme: ',usepingpong write(6,'(a24,1x,es8.1,/)')' unitlength: ',numerics_unitlength @@ -266,7 +263,6 @@ subroutine numerics_init write(6,'(a24,1x,es8.1)') ' err_curl_tolRel: ',err_curl_tolRel write(6,'(a24,1x,es8.1)') ' polarAlpha: ',polarAlpha write(6,'(a24,1x,es8.1)') ' polarBeta: ',polarBeta - write(6,'(a24,1x,a)') ' PETSc_options: ',trim(petsc_options) #endif !-------------------------------------------------------------------------------------------------- @@ -274,16 +270,17 @@ subroutine numerics_init #ifdef FEM write(6,'(a24,1x,i8)') ' integrationOrder: ',integrationOrder write(6,'(a24,1x,i8)') ' structOrder: ',structOrder - write(6,'(a24,1x,a)') ' PETSc_options: ',trim(petsc_options) write(6,'(a24,1x,L8)') ' B-Bar stabilisation: ',BBarStabilisation #endif +#ifdef PETSC + write(6,'(a24,1x,a)') ' PETSc_options: ',trim(petsc_options) +#endif + !-------------------------------------------------------------------------------------------------- ! sanity checks if (defgradTolerance <= 0.0_pReal) call IO_error(301,ext_msg='defgradTolerance') if (iJacoStiffness < 1) call IO_error(301,ext_msg='iJacoStiffness') - if (numerics_integrator <= 0 .or. numerics_integrator >= 6) & - call IO_error(301,ext_msg='integrator') if (numerics_unitlength <= 0.0_pReal) call IO_error(301,ext_msg='unitlength') if (residualStiffness < 0.0_pReal) call IO_error(301,ext_msg='residualStiffness') if (itmax <= 1) call IO_error(301,ext_msg='itmax') From 97e89f3f88ec4cb67e76160c5c51d1ac27e544fa Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 11:09:30 +0200 Subject: [PATCH 12/23] nonlocal can run in local mode --- src/constitutive_plastic_nonlocal.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index 35cc7bdd1..99ccf9ec4 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -320,6 +320,8 @@ module subroutine plastic_nonlocal_init prm%fEdgeMultiplication = config%getFloat('edgemultiplication') prm%shortRangeStressCorrection = config%keyExists('/shortrangestresscorrection/') + plasticState(p)%nonlocal = config%KeyExists('/nonlocal/') + !-------------------------------------------------------------------------------------------------- ! sanity checks if (any(prm%burgers < 0.0_pReal)) extmsg = trim(extmsg)//' burgers' @@ -386,7 +388,6 @@ module subroutine plastic_nonlocal_init call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,sizeDeltaState) - plasticState(p)%nonlocal = .true. plasticState(p)%offsetDeltaState = 0 ! ToDo: state structure does not follow convention st0%rho => plasticState(p)%state0 (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:) From 5af53f0be75617d1848f3135ddaf9ab534e922fa Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 11:30:01 +0200 Subject: [PATCH 13/23] nonlocal is a property of the phase --- src/crystallite.f90 | 99 ++++++++++++++++++--------------------------- 1 file changed, 40 insertions(+), 59 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 7b3ee4115..611a7551d 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1029,9 +1029,10 @@ subroutine integrateStateFPI(todo) 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) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then + p = material_phaseAt(g,e) + if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then - p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) + c = material_phaseMemberAt(g,i,e) broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & crystallite_partionedF0, & @@ -1039,8 +1040,7 @@ subroutine integrateStateFPI(todo) crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e) - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle sizeDotState = plasticState(p)%sizeDotState @@ -1112,8 +1112,7 @@ subroutine integrateStateFPI(todo) endif enddo iteration - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. endif enddo; enddo; enddo @@ -1169,9 +1168,10 @@ subroutine integrateStateEuler(todo) 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) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then + p = material_phaseAt(g,e) + if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then - p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) + c = material_phaseMemberAt(g,i,e) broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & crystallite_partionedF0, & @@ -1179,8 +1179,7 @@ subroutine integrateStateEuler(todo) crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e) - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle sizeDotState = plasticState(p)%sizeDotState @@ -1195,13 +1194,11 @@ subroutine integrateStateEuler(todo) enddo broken = stateJump(g,i,e) - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle broken = integrateStress(g,i,e) - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. crystallite_converged(g,i,e) = .not. broken @@ -1240,9 +1237,10 @@ subroutine integrateStateAdaptiveEuler(todo) 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) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then + p = material_phaseAt(g,e) + if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then - p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) + c = material_phaseMemberAt(g,i,e) broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & crystallite_partionedF0, & @@ -1250,8 +1248,7 @@ subroutine integrateStateAdaptiveEuler(todo) crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e) - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle sizeDotState = plasticState(p)%sizeDotState @@ -1269,13 +1266,11 @@ subroutine integrateStateAdaptiveEuler(todo) enddo broken = stateJump(g,i,e) - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle broken = integrateStress(g,i,e) - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & @@ -1284,8 +1279,7 @@ subroutine integrateStateAdaptiveEuler(todo) crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e) - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1322,16 +1316,16 @@ subroutine integrateStateRK4(todo) logical, dimension(:,:,:), intent(in) :: todo - real(pReal), dimension(3,3), parameter :: & - A = reshape([& + real(pReal), dimension(3,3), parameter :: & + A = reshape([& 0.5_pReal, 0.0_pReal, 0.0_pReal, & 0.0_pReal, 0.5_pReal, 0.0_pReal, & 0.0_pReal, 0.0_pReal, 1.0_pReal], & [3,3]) - real(pReal), dimension(3), parameter :: & - CC = [0.5_pReal, 0.5_pReal, 1.0_pReal] ! factor giving the fraction of the original timestep used for Runge Kutta Integration - 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] ! weight of slope used for Runge Kutta integration (final weight divided by 6) + real(pReal), dimension(3), parameter :: & + CC = [0.5_pReal, 0.5_pReal, 1.0_pReal] ! factor giving the fraction of the original timestep used for Runge Kutta Integration + 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] ! weight of slope used for Runge Kutta integration (final weight divided by 6) integer :: & e, & ! element index in element loop @@ -1353,9 +1347,11 @@ subroutine integrateStateRK4(todo) 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) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then + p = material_phaseAt(g,e) + if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then + + 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, & @@ -1363,8 +1359,7 @@ subroutine integrateStateRK4(todo) crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e) - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle do stage = 1,3 @@ -1400,8 +1395,6 @@ subroutine integrateStateRK4(todo) enddo broken = integrateStress(g,i,e,CC(stage)) - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. if(broken) exit broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & @@ -1409,13 +1402,10 @@ subroutine integrateStateRK4(todo) crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e)*CC(stage), g,i,e) - - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. if(broken) exit enddo - + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle sizeDotState = plasticState(p)%sizeDotState @@ -1439,13 +1429,11 @@ subroutine integrateStateRK4(todo) enddo broken = stateJump(g,i,e) - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle broken = integrateStress(g,i,e) - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. crystallite_converged(g,i,e) = .not. broken endif @@ -1505,9 +1493,10 @@ subroutine integrateStateRKCK45(todo) 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) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then + p = material_phaseAt(g,e) + if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then - p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) + c = material_phaseMemberAt(g,i,e) broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & crystallite_partionedF0, & @@ -1515,8 +1504,7 @@ subroutine integrateStateRKCK45(todo) crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e) - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle do stage = 1,5 @@ -1552,8 +1540,6 @@ subroutine integrateStateRKCK45(todo) enddo broken = integrateStress(g,i,e,CC(stage)) - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. if(broken) exit broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & @@ -1561,12 +1547,10 @@ subroutine integrateStateRKCK45(todo) crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e)*CC(stage), g,i,e) - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. if(broken) exit enddo - + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle sizeDotState = plasticState(p)%sizeDotState @@ -1595,18 +1579,15 @@ subroutine integrateStateRKCK45(todo) sourceState(p)%p(s)%state(1:sizeDotState,c), & sourceState(p)%p(s)%atol(1:sizeDotState)) enddo - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle broken = stateJump(g,i,e) - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle broken = integrateStress(g,i,e) - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. crystallite_converged(g,i,e) = .not. broken endif From b5efaa08a47ea7e0a662b7b7cbf3ca8e4762aa69 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 12:55:49 +0200 Subject: [PATCH 14/23] use already known mappings --- src/constitutive.f90 | 35 ++++++++++++++++++----------------- src/crystallite.f90 | 35 ++++++++++++++++------------------- 2 files changed, 34 insertions(+), 36 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 582222dae..e7f5628ed 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -709,12 +709,14 @@ end subroutine constitutive_hooke_SandItsTangents !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el) result(broken) +function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) result(broken) integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point - el !< element + el, & !< element + phase, & + of real(pReal), intent(in) :: & subdt !< timestep real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & @@ -727,17 +729,14 @@ function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el real(pReal), dimension(3,3) :: & Mp integer :: & - phase, & ho, & !< homogenization tme, & !< thermal member position i, & !< counter in source loop - instance, of + instance logical :: broken ho = material_homogenizationAt(el) tme = thermalMapping(ho)%p(ip,el) - of = material_phasememberAt(ipc,ip,el) - phase = material_phaseAt(ipc,el) instance = phase_plasticityInstance(phase) Mp = matmul(matmul(transpose(Fi),Fi),S) @@ -794,12 +793,14 @@ end function constitutive_collectDotState !> @brief for constitutive models having an instantaneous change of state !> will return false if delta state is not needed/supported by the constitutive model !-------------------------------------------------------------------------------------------------- -function constitutive_deltaState(S, Fe, Fi, ipc, ip, el) result(broken) +function constitutive_deltaState(S, Fe, Fi, ipc, ip, el, phase, of) result(broken) integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point - el !< element + el, & !< element + phase, & + of real(pReal), intent(in), dimension(3,3) :: & S, & !< 2nd Piola Kirchhoff stress Fe, & !< elastic deformation gradient @@ -808,27 +809,28 @@ function constitutive_deltaState(S, Fe, Fi, ipc, ip, el) result(broken) Mp integer :: & i, & - instance, of, & - phase + instance logical :: & broken - Mp = matmul(matmul(transpose(Fi),Fi),S) - of = material_phasememberAt(ipc,ip,el) - phase = material_phaseAt(ipc,el) - instance = phase_plasticityInstance(material_phaseAt(ipc,el)) + instance = phase_plasticityInstance(phase) plasticityType: select case (phase_plasticity(phase)) case (PLASTICITY_KINEHARDENING_ID) plasticityType call plastic_kinehardening_deltaState(Mp,instance,of) + broken = any(IEEE_is_NaN(plasticState(phase)%deltaState(:,of))) case (PLASTICITY_NONLOCAL_ID) plasticityType call plastic_nonlocal_deltaState(Mp,instance,of,ip,el) + broken = any(IEEE_is_NaN(plasticState(phase)%deltaState(:,of))) + + case default + broken = .false. end select plasticityType - broken = any(IEEE_is_NaN(plasticState(phase)%deltaState(:,of))) + sourceLoop: do i = 1, phase_Nsources(phase) @@ -837,11 +839,10 @@ function constitutive_deltaState(S, Fe, Fi, ipc, ip, el) result(broken) case (SOURCE_damage_isoBrittle_ID) sourceType call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(ipc,ip,el), Fe, & ipc, ip, el) + broken = broken .or. any(IEEE_is_NaN(sourceState(phase)%p(i)%deltaState(:,of))) end select sourceType - broken = broken .or. any(IEEE_is_NaN(sourceState(phase)%p(i)%deltaState(:,of))) - enddo SourceLoop end function constitutive_deltaState diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 611a7551d..054bc2105 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1038,7 +1038,7 @@ subroutine integrateStateFPI(todo) crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e) + crystallite_subdt(g,i,e), g,i,e,p,c) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1073,7 +1073,7 @@ subroutine integrateStateFPI(todo) crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e) + crystallite_subdt(g,i,e), g,i,e,p,c) if(broken) exit iteration sizeDotState = plasticState(p)%sizeDotState @@ -1107,7 +1107,7 @@ subroutine integrateStateFPI(todo) enddo if(crystallite_converged(g,i,e)) then - broken = stateJump(g,i,e) + broken = stateJump(g,i,e,p,c) exit iteration endif @@ -1177,7 +1177,7 @@ subroutine integrateStateEuler(todo) crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e) + crystallite_subdt(g,i,e), g,i,e,p,c) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1193,7 +1193,7 @@ subroutine integrateStateEuler(todo) * crystallite_subdt(g,i,e) enddo - broken = stateJump(g,i,e) + broken = stateJump(g,i,e,p,c) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1246,7 +1246,7 @@ subroutine integrateStateAdaptiveEuler(todo) crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e) + crystallite_subdt(g,i,e), g,i,e,p,c) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1265,7 +1265,7 @@ subroutine integrateStateAdaptiveEuler(todo) + sourceState(p)%p(s)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) enddo - broken = stateJump(g,i,e) + broken = stateJump(g,i,e,p,c) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1277,7 +1277,7 @@ subroutine integrateStateAdaptiveEuler(todo) crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e) + crystallite_subdt(g,i,e), g,i,e,p,c) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1357,7 +1357,7 @@ subroutine integrateStateRK4(todo) crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e) + crystallite_subdt(g,i,e), g,i,e,p,c) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1401,7 +1401,7 @@ subroutine integrateStateRK4(todo) crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & - crystallite_subdt(g,i,e)*CC(stage), g,i,e) + crystallite_subdt(g,i,e)*CC(stage), g,i,e,p,c) if(broken) exit enddo @@ -1428,7 +1428,7 @@ subroutine integrateStateRK4(todo) * crystallite_subdt(g,i,e) enddo - broken = stateJump(g,i,e) + broken = stateJump(g,i,e,p,c) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1502,7 +1502,7 @@ subroutine integrateStateRKCK45(todo) crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e) + crystallite_subdt(g,i,e), g,i,e,p,c) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1546,7 +1546,7 @@ subroutine integrateStateRKCK45(todo) crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & - crystallite_subdt(g,i,e)*CC(stage), g,i,e) + crystallite_subdt(g,i,e)*CC(stage), g,i,e,p,c) if(broken) exit enddo @@ -1582,7 +1582,7 @@ subroutine integrateStateRKCK45(todo) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle - broken = stateJump(g,i,e) + broken = stateJump(g,i,e,p,c) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1631,7 +1631,7 @@ end function converged !> @brief calculates a jump in the state according to the current state and the current stress !> returns true, if state jump was successfull or not needed. false indicates NaN in delta state !-------------------------------------------------------------------------------------------------- -function stateJump(ipc,ip,el) result(broken) +function stateJump(ipc,ip,el,p,c) result(broken) integer, intent(in):: & el, & ! element index @@ -1646,13 +1646,10 @@ function stateJump(ipc,ip,el) result(broken) mySize logical :: broken - c = material_phaseMemberAt(ipc,ip,el) - p = material_phaseAt(ipc,el) - broken = constitutive_deltaState(crystallite_S(1:3,1:3,ipc,ip,el), & crystallite_Fe(1:3,1:3,ipc,ip,el), & crystallite_Fi(1:3,1:3,ipc,ip,el), & - ipc,ip,el) + ipc,ip,el,p,c) if(broken) return myOffset = plasticState(p)%offsetDeltaState From 70dd06c4ecca2e9df9fe3ce19c7618dd0cc1aacc Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 14:42:38 +0200 Subject: [PATCH 15/23] constitutive should handle state jump alone --- src/constitutive.f90 | 21 ++++++++++++++- src/crystallite.f90 | 62 +++++++++++--------------------------------- 2 files changed, 35 insertions(+), 48 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index e7f5628ed..e2c9dbc05 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -809,7 +809,9 @@ function constitutive_deltaState(S, Fe, Fi, ipc, ip, el, phase, of) result(broke Mp integer :: & i, & - instance + instance, & + myOffset, & + mySize logical :: & broken @@ -831,6 +833,17 @@ function constitutive_deltaState(S, Fe, Fi, ipc, ip, el, phase, of) result(broke end select plasticityType + if(.not. broken) then + select case(phase_plasticity(phase)) + case (PLASTICITY_NONLOCAL_ID,PLASTICITY_KINEHARDENING_ID) + + myOffset = plasticState(phase)%offsetDeltaState + mySize = plasticState(phase)%sizeDeltaState + plasticState(phase)%state(myOffset + 1:myOffset + mySize,of) = & + plasticState(phase)%state(myOffset + 1:myOffset + mySize,of) + plasticState(phase)%deltaState(1:mySize,of) + end select + endif + sourceLoop: do i = 1, phase_Nsources(phase) @@ -840,6 +853,12 @@ function constitutive_deltaState(S, Fe, Fi, ipc, ip, el, phase, of) result(broke call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(ipc,ip,el), Fe, & ipc, ip, el) broken = broken .or. any(IEEE_is_NaN(sourceState(phase)%p(i)%deltaState(:,of))) + if(.not. broken) then + myOffset = sourceState(phase)%p(i)%offsetDeltaState + mySize = sourceState(phase)%p(i)%sizeDeltaState + sourceState(phase)%p(i)%state(myOffset + 1: myOffset + mySize,of) = & + sourceState(phase)%p(i)%state(myOffset + 1: myOffset + mySize,of) + sourceState(phase)%p(i)%deltaState(1:mySize,of) + endif end select sourceType diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 054bc2105..8bceff9a4 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1107,7 +1107,9 @@ subroutine integrateStateFPI(todo) enddo if(crystallite_converged(g,i,e)) then - broken = stateJump(g,i,e,p,c) + 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 @@ -1193,7 +1195,9 @@ subroutine integrateStateEuler(todo) * crystallite_subdt(g,i,e) enddo - broken = stateJump(g,i,e,p,c) + 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 @@ -1265,7 +1269,9 @@ subroutine integrateStateAdaptiveEuler(todo) + sourceState(p)%p(s)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) enddo - broken = stateJump(g,i,e,p,c) + 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 @@ -1428,7 +1434,9 @@ subroutine integrateStateRK4(todo) * crystallite_subdt(g,i,e) enddo - broken = stateJump(g,i,e,p,c) + 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 @@ -1582,7 +1590,9 @@ subroutine integrateStateRKCK45(todo) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle - broken = stateJump(g,i,e,p,c) + 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 @@ -1627,48 +1637,6 @@ logical pure function converged(residuum,state,atol) end function converged -!-------------------------------------------------------------------------------------------------- -!> @brief calculates a jump in the state according to the current state and the current stress -!> returns true, if state jump was successfull or not needed. false indicates NaN in delta state -!-------------------------------------------------------------------------------------------------- -function stateJump(ipc,ip,el,p,c) result(broken) - - integer, intent(in):: & - el, & ! element index - ip, & ! integration point index - ipc ! grain index - - integer :: & - c, & - p, & - mySource, & - myOffset, & - mySize - logical :: broken - - broken = constitutive_deltaState(crystallite_S(1:3,1:3,ipc,ip,el), & - crystallite_Fe(1:3,1:3,ipc,ip,el), & - crystallite_Fi(1:3,1:3,ipc,ip,el), & - ipc,ip,el,p,c) - if(broken) return - - myOffset = plasticState(p)%offsetDeltaState - mySize = plasticState(p)%sizeDeltaState - - - plasticState(p)%state(myOffset + 1:myOffset + mySize,c) = & - plasticState(p)%state(myOffset + 1:myOffset + mySize,c) + plasticState(p)%deltaState(1:mySize,c) - - do mySource = 1, phase_Nsources(p) - myOffset = sourceState(p)%p(mySource)%offsetDeltaState - mySize = sourceState(p)%p(mySource)%sizeDeltaState - sourceState(p)%p(mySource)%state(myOffset + 1: myOffset + mySize,c) = & - sourceState(p)%p(mySource)%state(myOffset + 1: myOffset + mySize,c) + sourceState(p)%p(mySource)%deltaState(1:mySize,c) - enddo - -end function stateJump - - !-------------------------------------------------------------------------------------------------- !> @brief Write current restart information (Field and constitutive data) to file. ! ToDo: Merge data into one file for MPI, move state to constitutive and homogenization, respectively From cf5fcf389b1502da09e86826684563c5391e4d40 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 17:53:07 +0200 Subject: [PATCH 16/23] phase is a property of the element and we have no homogenization for nonlocal --- src/crystallite.f90 | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 8bceff9a4..967a34dbb 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -615,14 +615,16 @@ subroutine crystallite_orientations enddo; enddo; enddo !$OMP END PARALLEL DO - nonlocalPresent: if (any(plasticState%nonLocal)) then + nonlocalPresent: if (any(plasticState%nonlocal)) then !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1),FEsolving_execIP(2) - if (plasticState(material_phaseAt(1,e))%nonLocal) & + if (plasticState(material_phaseAt(1,e))%nonlocal) then + do i = FEsolving_execIP(1),FEsolving_execIP(2) call plastic_nonlocal_updateCompatibility(crystallite_orientation, & phase_plasticityInstance(material_phaseAt(i,e)),i,e) - enddo; enddo + enddo + endif + enddo !$OMP END PARALLEL DO endif nonlocalPresent From d0d963a2cc8f368aca29ee2dcfd610600331051d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 17:58:50 +0200 Subject: [PATCH 17/23] set independent of number of slip systems --- src/constitutive_plastic_nonlocal.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index 99ccf9ec4..d530f6651 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -320,7 +320,6 @@ module subroutine plastic_nonlocal_init prm%fEdgeMultiplication = config%getFloat('edgemultiplication') prm%shortRangeStressCorrection = config%keyExists('/shortrangestresscorrection/') - plasticState(p)%nonlocal = config%KeyExists('/nonlocal/') !-------------------------------------------------------------------------------------------------- ! sanity checks @@ -388,6 +387,7 @@ module subroutine plastic_nonlocal_init call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,sizeDeltaState) + plasticState(p)%nonlocal = config%KeyExists('/nonlocal/') plasticState(p)%offsetDeltaState = 0 ! ToDo: state structure does not follow convention st0%rho => plasticState(p)%state0 (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:) From 3a4bb59057bcaeeac4f476a7ef992f64f0825bea Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 18:58:48 +0200 Subject: [PATCH 18/23] no need to store the same information multiple times --- src/crystallite.f90 | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 967a34dbb..0d905aaa2 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -69,8 +69,7 @@ module crystallite logical, dimension(:,:,:), allocatable, public :: & crystallite_requested !< used by upper level (homogenization) to request crystallite calculation logical, dimension(:,:,:), allocatable :: & - crystallite_converged, & !< convergence flag - crystallite_localPlasticity !< indicates this grain to have purely local constitutive law + crystallite_converged !< convergence flag type :: tOutput !< new requested output (per phase) character(len=pStringLen), allocatable, dimension(:) :: & @@ -158,7 +157,6 @@ subroutine crystallite_init allocate(crystallite_orientation(cMax,iMax,eMax)) - allocate(crystallite_localPlasticity(cMax,iMax,eMax), source=.true.) allocate(crystallite_requested(cMax,iMax,eMax), source=.false.) allocate(crystallite_converged(cMax,iMax,eMax), source=.true.) @@ -238,7 +236,6 @@ subroutine crystallite_init / math_det33(crystallite_Fp0(1:3,1:3,c,i,e))**(1.0_pReal/3.0_pReal) crystallite_Fi0(1:3,1:3,c,i,e) = constitutive_initialFi(c,i,e) crystallite_F0(1:3,1:3,c,i,e) = math_I3 - crystallite_localPlasticity(c,i,e) = phase_localPlasticity(material_phaseAt(c,e)) crystallite_Fe(1:3,1:3,c,i,e) = math_inv33(matmul(crystallite_Fi0(1:3,1:3,c,i,e), & crystallite_Fp0(1:3,1:3,c,i,e))) ! assuming that euler angles are given in internal strain free configuration crystallite_Fp(1:3,1:3,c,i,e) = crystallite_Fp0(1:3,1:3,c,i,e) @@ -248,7 +245,7 @@ subroutine crystallite_init enddo !$OMP END PARALLEL DO - if(any(.not. crystallite_localPlasticity) .and. .not. usePingPong) call IO_error(601) ! exit if nonlocal but no ping-pong ToDo: Why not check earlier? or in nonlocal? + if(any(plasticState%nonlocal) .and. .not. usePingPong) call IO_error(601) ! exit if nonlocal but no ping-pong ToDo: Why not check earlier? or in nonlocal? crystallite_partionedFp0 = crystallite_Fp0 crystallite_partionedFi0 = crystallite_Fi0 @@ -277,7 +274,6 @@ subroutine crystallite_init write(6,'(a42,1x,i10)') ' # of elements: ', eMax write(6,'(a42,1x,i10)') ' # of integration points/element: ', iMax write(6,'(a42,1x,i10)') 'max # of constituents/integration point: ', cMax - write(6,'(a42,1x,i10)') ' # of nonlocal constituents: ',count(.not. crystallite_localPlasticity) flush(6) endif @@ -1606,7 +1602,7 @@ subroutine integrateStateRKCK45(todo) enddo; enddo; enddo !$OMP END PARALLEL DO - if (nonlocalBroken) call nonlocalConvergenceCheck + if(nonlocalBroken) call nonlocalConvergenceCheck end subroutine integrateStateRKCK45 @@ -1617,7 +1613,14 @@ end subroutine integrateStateRKCK45 !-------------------------------------------------------------------------------------------------- subroutine nonlocalConvergenceCheck - where( .not. crystallite_localPlasticity) crystallite_converged = .false. + integer :: e,i,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) crystallite_converged(1,i,e) = .false. + enddo + enddo end subroutine nonlocalConvergenceCheck From d9806cb7f3ffe9a031b7946f3fbb33c193d4db26 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 2 Apr 2020 00:23:59 +0200 Subject: [PATCH 19/23] do not clutter with nonlocal checks --- src/crystallite.f90 | 30 ++++++------------------------ 1 file changed, 6 insertions(+), 24 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 0d905aaa2..7fc51b1c9 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1037,7 +1037,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 @@ -1113,7 +1112,6 @@ subroutine integrateStateFPI(todo) enddo iteration if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. - endif enddo; enddo; enddo !$OMP END PARALLEL DO @@ -1178,7 +1176,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 @@ -1201,9 +1198,7 @@ subroutine integrateStateEuler(todo) broken = integrateStress(g,i,e) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. - crystallite_converged(g,i,e) = .not. broken - endif enddo; enddo; enddo !$OMP END PARALLEL DO @@ -1239,6 +1234,7 @@ subroutine integrateStateAdaptiveEuler(todo) 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 @@ -1249,8 +1245,6 @@ 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 .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle sizeDotState = plasticState(p)%sizeDotState @@ -1270,11 +1264,9 @@ 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 .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle broken = integrateStress(g,i,e) - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & @@ -1282,13 +1274,10 @@ 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 .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle 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), & @@ -1296,7 +1285,6 @@ subroutine integrateStateAdaptiveEuler(todo) 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), & @@ -1305,6 +1293,7 @@ subroutine integrateStateAdaptiveEuler(todo) enddo endif + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. enddo; enddo; enddo !$OMP END PARALLEL DO @@ -1351,6 +1340,7 @@ subroutine integrateStateRK4(todo) 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 @@ -1362,8 +1352,6 @@ subroutine integrateStateRK4(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 do stage = 1,3 @@ -1409,7 +1397,6 @@ subroutine integrateStateRK4(todo) if(broken) exit enddo - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle sizeDotState = plasticState(p)%sizeDotState @@ -1435,14 +1422,13 @@ subroutine integrateStateRK4(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 + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. enddo; enddo; enddo !$OMP END PARALLEL DO @@ -1499,6 +1485,7 @@ subroutine integrateStateRKCK45(todo) 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 @@ -1509,8 +1496,6 @@ subroutine integrateStateRKCK45(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 do stage = 1,5 @@ -1556,7 +1541,6 @@ subroutine integrateStateRKCK45(todo) if(broken) exit enddo - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle sizeDotState = plasticState(p)%sizeDotState @@ -1585,20 +1569,18 @@ subroutine integrateStateRKCK45(todo) sourceState(p)%p(s)%state(1:sizeDotState,c), & sourceState(p)%p(s)%atol(1:sizeDotState)) enddo - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle 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 + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. enddo; enddo; enddo !$OMP END PARALLEL DO From 565cf8239f077c8830dd195cf41bef3ae99fcd51 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 2 Apr 2020 09:29:58 +0200 Subject: [PATCH 20/23] can be done in parallel --- src/crystallite.f90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 7fc51b1c9..595947803 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1597,12 +1597,14 @@ subroutine nonlocalConvergenceCheck integer :: e,i,p + !$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) crystallite_converged(1,i,e) = .false. enddo enddo + !$OMP END PARALLEL DO end subroutine nonlocalConvergenceCheck From b375af83a3dcedc8dae6ef564ab0c1800d396a7a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 10 Apr 2020 15:52:54 +0200 Subject: [PATCH 21/23] bugfix for issue introduced in last merge --- src/crystallite.f90 | 55 +-------------------------------------------- 1 file changed, 1 insertion(+), 54 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 743bcf73b..f0b5ecf08 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -301,7 +301,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) e, & !< counter in element loop startIP, endIP, & s - logical, dimension(homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: todo !ToDo: need to set some values to false in has of different Ngrains + logical, dimension(homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: todo !ToDo: need to set some values to false for different Ngrains #ifdef DEBUG if (iand(debug_level(debug_crystallite),debug_levelSelective) /= 0 & @@ -1626,59 +1626,6 @@ logical pure function converged(residuum,state,atol) end function converged -!-------------------------------------------------------------------------------------------------- -!> @brief calculates a jump in the state according to the current state and the current stress -!> returns true, if state jump was successfull or not needed. false indicates NaN in delta state -!-------------------------------------------------------------------------------------------------- -logical function stateJump(ipc,ip,el) - - integer, intent(in):: & - el, & ! element index - ip, & ! integration point index - ipc ! grain index - - integer :: & - c, & - p, & - mySource, & - myOffset, & - mySize - - c = material_phaseMemberAt(ipc,ip,el) - p = material_phaseAt(ipc,el) - - call constitutive_collectDeltaState(crystallite_S(1:3,1:3,ipc,ip,el), & - crystallite_Fe(1:3,1:3,ipc,ip,el), & - crystallite_Fi(1:3,1:3,ipc,ip,el), & - ipc,ip,el) - - myOffset = plasticState(p)%offsetDeltaState - mySize = plasticState(p)%sizeDeltaState - - if( any(IEEE_is_NaN(plasticState(p)%deltaState(1:mySize,c)))) then - stateJump = .false. - return - endif - - plasticState(p)%state(myOffset + 1:myOffset + mySize,c) = & - plasticState(p)%state(myOffset + 1:myOffset + mySize,c) + plasticState(p)%deltaState(1:mySize,c) - - do mySource = 1, phase_Nsources(p) - myOffset = sourceState(p)%p(mySource)%offsetDeltaState - mySize = sourceState(p)%p(mySource)%sizeDeltaState - if (any(IEEE_is_NaN(sourceState(p)%p(mySource)%deltaState(1:mySize,c)))) then - stateJump = .false. - return - endif - sourceState(p)%p(mySource)%state(myOffset + 1: myOffset + mySize,c) = & - sourceState(p)%p(mySource)%state(myOffset + 1: myOffset + mySize,c) + sourceState(p)%p(mySource)%deltaState(1:mySize,c) - enddo - - stateJump = .true. - -end function stateJump - - !-------------------------------------------------------------------------------------------------- !> @brief Write current restart information (Field and constitutive data) to file. ! ToDo: Merge data into one file for MPI, move state to constitutive and homogenization, respectively From 399a0f1b66cd771b60f7f61c113227ea9018e255 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 10 Apr 2020 18:53:59 +0200 Subject: [PATCH 22/23] store information about state size only once --- src/crystallite.f90 | 79 ++++++++++++++++++++++----------------------- 1 file changed, 39 insertions(+), 40 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index f0b5ecf08..f405fc930 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1012,7 +1012,9 @@ subroutine integrateStateFPI(todo) p, & c, & s, & - sizeDotState + size_pl + integer, dimension(maxval(phase_Nsources)) :: & + size_so real(pReal) :: & zeta real(pReal), dimension(max(constitutive_plasticity_maxSizeDotState,constitutive_source_maxSizeDotState)) :: & @@ -1023,7 +1025,7 @@ subroutine integrateStateFPI(todo) nonlocalBroken, broken nonlocalBroken = .false. - !$OMP PARALLEL DO PRIVATE(sizeDotState,r,zeta,p,c,plastic_dotState_p1, plastic_dotState_p2,source_dotState,broken) + !$OMP PARALLEL DO PRIVATE(size_pl,size_so,r,zeta,p,c,plastic_dotState_p1, plastic_dotState_p2,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)) @@ -1040,17 +1042,17 @@ subroutine integrateStateFPI(todo) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle - 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) - plastic_dotState_p2 = 0.0_pReal * plasticState(p)%dotState (1:sizeDotState,c) ! ToDo can be done smarter/clearer + 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_p2 = 0.0_pReal * plasticState(p)%dotState (1:size_pl,c) 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) - source_dotState(1:sizeDotState,2,s) = 0.0_pReal + 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 @@ -1058,49 +1060,46 @@ subroutine integrateStateFPI(todo) if(nIterationState > 1) plastic_dotState_p2 = plastic_dotState_p1 plastic_dotState_p1 = plasticState(p)%dotState(:,c) do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState - if(nIterationState > 1) source_dotState(1:sizeDotState,2,s) = source_dotState(1:sizeDotState,1,s) - source_dotState(1:sizeDotState,1,s) = sourceState(p)%p(s)%dotState(:,c) + 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 = 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) + 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 - sizeDotState = plasticState(p)%sizeDotState zeta = damper(plasticState(p)%dotState(:,c),plastic_dotState_p1,plastic_dotState_p2) plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * zeta & + plastic_dotState_p1 * (1.0_pReal - zeta) - r(1:SizeDotState) = plasticState(p)%state (1:sizeDotState,c) & - - plasticState(p)%subState0(1:sizeDotState,c) & - - plasticState(p)%dotState (1:sizeDotState,c) * crystallite_subdt(g,i,e) - plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%state(1:sizeDotState,c) & - - r(1:sizeDotState) - crystallite_converged(g,i,e) = converged(r(1:sizeDotState), & - plasticState(p)%state(1:sizeDotState,c), & - plasticState(p)%atol(1:sizeDotState)) + 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) - sizeDotState = sourceState(p)%p(s)%sizeDotState zeta = damper(sourceState(p)%p(s)%dotState(:,c), & - source_dotState(1:sizeDotState,1,s),& - source_dotState(1:sizeDotState,2,s)) + 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:sizeDotState,1,s)* (1.0_pReal - zeta) - r(1: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) - sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%state(1:sizeDotState,c) & - - r(1:sizeDotState) + + 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:sizeDotState), & - sourceState(p)%p(s)%state(1:sizeDotState,c), & - sourceState(p)%p(s)%atol(1:sizeDotState)) + 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 From 1d5b1a17cd2795de3cb80eaafbbc6bafaebd0020 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 10 Apr 2020 19:36:29 +0200 Subject: [PATCH 23/23] same style for plastic and source --- src/crystallite.f90 | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index f405fc930..2d4b9785d 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1019,13 +1019,14 @@ subroutine integrateStateFPI(todo) zeta real(pReal), dimension(max(constitutive_plasticity_maxSizeDotState,constitutive_source_maxSizeDotState)) :: & r ! state residuum - real(pReal), dimension(:), allocatable :: plastic_dotState_p1, plastic_dotState_p2 + real(pReal), dimension(constitutive_plasticity_maxSizeDotState,2) :: & + plastic_dotState real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState logical :: & nonlocalBroken, broken nonlocalBroken = .false. - !$OMP PARALLEL DO PRIVATE(size_pl,size_so,r,zeta,p,c,plastic_dotState_p1, plastic_dotState_p2,source_dotState,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)) @@ -1046,7 +1047,7 @@ subroutine integrateStateFPI(todo) 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_p2 = 0.0_pReal * plasticState(p)%dotState (1:size_pl,c) + 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) & @@ -1057,8 +1058,8 @@ subroutine integrateStateFPI(todo) iteration: do NiterationState = 1, num%nState - if(nIterationState > 1) plastic_dotState_p2 = plastic_dotState_p1 - plastic_dotState_p1 = plasticState(p)%dotState(:,c) + 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) @@ -1074,9 +1075,10 @@ subroutine integrateStateFPI(todo) crystallite_subdt(g,i,e), g,i,e,p,c) if(broken) exit iteration - zeta = damper(plasticState(p)%dotState(:,c),plastic_dotState_p1,plastic_dotState_p2) + 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_p1 * (1.0_pReal - 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)